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This  is  an  interdisciplinary  project.  The  main  results  of  the  project  are: 


1. The  analytical  proof  of  the  global  convergence  property  using  a  sophisticated  mathematical  apparatus. 

2.  The  development  of  a  sophisticated  analytical  apparatus  for  establishing  the  relaxation  property  of  the  adaptivity  technique. 
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7.  Verification  of  the  globally  convergent  numerical  method  on  backscattering  experimental  data  for  targets  buried  in  the  ground.  This  case  is 
much  more  complicated  than  the  case  of  targets  in  air. 
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Program  Manager  of  the  Numerical  Analysis  Program  of  ARO. 
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15.  The  first  solution  of  a  long  standing  problem  about  uniqueness  of  a  phaseless  3-d  inverse  problem  of  quantum  scattering.  This  was  an 
open  question  since  the  publication  in  1977  of  the  well  known  book  of  K.  Chadan  and  P.  Sabatier  "Inverse  Problems  in  Quantum  Scattering 
Theory",  Springer,  New  York. 
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The  target  application  of  this  project  is  in  the  standoff  detection  and  identification  of  explosives,  such  as  lEDs  and  antipersonnel  land  mines, 
using  electric  signals  propagations.  The  numerical  methods  developed  in  this  project  use  experimental  data  to  image  both  dielectric  constants 
and  shapes  of  small  inclusions  buried  in  a  sand  box.  Inclusions  mimic  explosives.  A  systematic  study  of  twenty  five  cases  was  performed. 

Conventional  least  squares  cost  functionals  for  solving  Coefficient  Inverse  Problems  (CIPs)  face  the  fundamental  obstacle:  the  phenomenon 
of  multiple  local  minima  and  ravines.  This  results,  in  turn  in  the  local  convergence  of  conventional  numerical  methods  for  CIPs.  Therefore, 
the  following  question  is  both  the  most  crucial  and  the  most  challenging  one  in  a  numerical  treatment  of  any  CIP:  How  to  rigorously  obtain  a 
good  approximation  for  the  unknown  coefficient  without  any  advanced  knowledge  of  a  small  neighborhood  of  the  solution?  A  numerical 
method  providing  such  an  approximation  is  called  globally  convergent.  As  soon  as  this  approximation  is  found,  some  conventional  locally 
convergent  numerical  methods  can  be  used  as  refinement  tools. 

This  project  has  been  focused  on  the  development  of  a  globally  convergent  numerical  method  for  a  hyperbolic  CIP.  The  crucial  advantage  of 
this  method  is  that  it  does  not  need  a  good  first  guess  for  the  solution,  unlike  all  existing  ones.  Both  analytical  and  numerical  issues  were 
addressed.  It  was  observed  that  the  best  is  a  two-stage  numerical  procedure:  On  the  first  stage  the  globally  convergent  numerical  method 
delivers  a  good  approximation  for  the  solution.  And  on  the  second  stage  the  adaptivity  technique  refines  the  solution  obtained  on  the  first 
stage. 
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Quantum  Scattering  Theory",  Springer,  New  York,  1977. 

16.  Publication  of  a  book  summarizing  some  results  of  the  project  up  to  2012.  This  book  was  published  by  Springer,  the  World 
most  prestigious  publisher  of  the  scientific  literature. 

1 /.Publication  of  two  survey  papers. 
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1  .Addressing  a  need  of  the  Army  via  successful  work  with  experimental  data  collected  by  the  Forward  Looking  Radar  of  US 
Army  Research  Laboratory  (ARL).  The  globally  convergent  method  of  this  project  was  used. 

2.Transfer  of  a  ready-to-use  software  to  ARL  resulting  from  item  #1 .  This  software  works  with  the  real  data  of  the  Forward 
Looking  Radar  of  ARL. 

S.The  use  of  experimental  data  of  Item  #1  for  a  comparison  of  performances  of  the  globally  convergent  numerical  method  of  this 
project  and  the  classical  Krein  equation  method.  It  was  established  that  while  the  first  method  works  well,  the  second  one  fails 
for  these  data. 

4.  Four  presentations  of  results  were  given  to  Mr.  Brian  Burns,  Drs.  Anders  Sullivan  and  Lam  Nguyen,  ARL  engineers,  and  to 
Dr.  Joseph  D.  Myers,  the  Program  Manager  of  the  Numerical  Analysis  Program  of  ARO. 

5. Three  joint  publications  with  ARL  engineers  Drs.  Anders  Sullivan  and  Lam  Nguyen. 
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51, 2937-2948,  2013. 
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Abstract 

This  is  an  interdisciplinary  project.  The  main  results  of  the  project  are: 

1.  The  analytical  proof  of  the  global  convergence  property  using  a  sophisticated  mathe¬ 
matical  apparatus. 

2.  The  development  of  a  sophisticated  analytical  apparatus  for  establishing  the  relaxation 
property  of  the  adaptivity  technique. 

3.  Numerical  implementations  of  resulting  algorithms. 

4.  Numerical  verihcations  of  resulting  algorithms  on  computationally  simulated  data. 

5.  Assembling  an  experimental  apparatus  in  Microwave  Laboratory  of  University  of  North 
Carolina  at  Charlotte. 

6.  Verihcation  of  the  globally  convergent  numerical  method  on  backscattering  experimen¬ 
tal  data  for  targets  standing  in  air.  Targets  mimic  explosives. 

7.  Verihcation  of  the  globally  convergent  numerical  method  on  backscattering  experimen¬ 
tal  data  for  targets  buried  in  the  ground.  This  case  is  much  more  complicated  than 
the  case  of  targets  in  air. 

8.  An  experimental  and  numerical  reconstruction  evidence  of  the  super  resolution  phe¬ 
nomenon. 

9.  Addressing  a  need  of  the  Army  via  successful  work  with  experimental  data  collected 
by  the  Forward  Looking  Radar  of  US  Army  Research  Laboratory  (ARL).  The  globally 
convergent  method  of  this  project  was  used. 

10.  Transfer  of  a  ready-to-use  software  to  ARL  resulting  from  item  ^^9.  This  software 
works  with  the  real  data  of  the  Forward  Looking  Radar  of  ARL. 

11.  The  use  of  experimental  data  of  item  ^^9  for  a  comparison  of  performances  of  the 
globally  convergent  numerical  method  of  this  project  and  the  classical  Krein  equation 
method.  It  was  established  that  while  the  first  method  works  well,  the  second  one  fails 
for  these  data. 
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12.  Four  (4)  presentations  to  Mr.  Brian  Burns,  Drs.  Anders  Sullivan  and  Lam  Nguyen, 
ARL  engineers,  and  to  Dr.  Joseph  D.  Myers,  the  Program  Manager  of  the  Numerical 
Analysis  Program  of  ARO. 

13.  Three  joint  publications  with  ARL  engineers  Drs.  Anders  Sullivan  and  Lam  Nguyen. 

14.  A  new  globally  convergent  numerical  method  based  on  the  Carleman  Weight  Function. 

15.  The  hrst  solution  of  a  long  standing  problem  about  uniqueness  of  a  phaseless  3-d  inverse 
problem  of  quantum  scattering.  This  was  an  open  question  since  the  publication  in 
1977  of  the  well  known  book  of  K.  Chadan  and  P.  Sabatier  ’’Inverse  Problems  in 
Quantum  Scattering  Theory” ,  Springer,  New  York. 

16.  Publication  of  a  book  summarizing  some  results  of  the  project.  This  book  was  pub¬ 
lished  by  Springer,  the  World  most  prestigious  publisher  of  the  scientific  literature. 

17.  Publication  of  twenty  two  papers.  Two  of  them  are  survey  papers. 

The  target  application  of  this  project  is  in  the  standoff  detection  and  identihcation  of 
explosives,  such  as  lEDs  and  antipersonnel  land  mines,  using  electric  signals  propagations. 
The  nnmerical  methods  developed  in  this  project  use  experimental  data  to  image  both  di¬ 
electric  constants  and  shapes  of  small  inclusions  buried  in  a  sand  box.  Inclusions  mimic 
explosives.  A  systematic  study  of  twenty  five  cases  was  performed. 

Conventional  least  sqnares  cost  fnnctionals  for  solving  Coefficient  Inverse  Problems  (CIPs) 
face  the  fnndamental  obstacle:  the  phenomenon  of  mnltiple  local  minima  and  ravines.  This 
results,  in  tnrn  in  the  local  convergence  of  conventional  nnmerical  methods  for  CIPs.  There¬ 
fore,  the  following  qnestion  is  both  the  most  crucial  and  the  most  challenging  one  in  a 
nnmerical  treatment  of  any  CIP:  How  to  rigorously  obtain  a  good  approximation  for  the  un¬ 
known  coefficient  without  any  advanced  knowledge  of  a  small  neighborhood  of  the  solution? 
A  nnmerical  method  providing  snch  an  approximation  is  called  globally  convergent.  As  soon 
as  this  approximation  is  fonnd,  some  conventional  locally  convergent  nnmerical  methods  can 
be  used  as  rehnement  tools. 

This  project  has  been  focused  on  the  development  of  a  globally  convergent  numerical 
method  for  a  hyperbolic  CIP.  The  crucial  advantage  of  this  method  is  that  it  does  not  need 
a  good  hrst  guess  for  the  solution,  unlike  all  existing  ones.  Both  analytical  and  numerical 
issues  were  addressed.  It  was  observed  that  the  best  is  a  two-stage  numerical  procedure:  On 
the  hrst  stage  the  globally  convergent  nnmerical  method  delivers  a  good  approximation  for 
the  solution.  And  on  the  second  stage  the  adaptivity  techniqne  rehnes  the  solntion  obtained 
on  the  hrst  stage. 
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1  Introduction 


This  is  an  interdisciplinary  pioiect.  It  combines  analytical  studies,  studies  of  computationally 
simulated  data  and  studies  of  experimental  data. 

The  target  application  of  this  project  is  detection  of  explosives,  such  as  lEDs,  plastic 
antipersonnel  land  mines,  etc.  using  electric  signals  propagations.  Indeed,  it  is  well  known 
from  tables  of  dielectric  constants  [65,  66]  that  dielectric  constants  in  explosives  are  higher 
than  those  of  regular  materials.  Hence,  the  idea  is  to  image  both  dielectric  constants  and 
shapes  of  small  sharp  inclusions  embedded  in  an  otherwise  slowly  changing  background. 

Different  imaging  methods  have  been  applied  to  this  type  of  measurements  to  obtain 
geometrical  information  such  as  shapes,  sizes  and  locations  of  the  targets,  see,  e.g.  [47]. 
However,  the  dielectric  constants,  which  characterize  the  targets  in  terms  of  their  constituent 
materials,  are  much  more  difficult  to  estimate.  Unlike  this,  the  globally  convergent  numerical 
method  of  the  current  project  can  image  both  dielectric  constants  and  shapes  of  inclusions 
mimicking  explosives.  The  crucial  advantage  of  this  method  is  that  it  does  not  need  a  good 
hrst  guess  for  the  solution,  unlike  all  existing  ones. 

A  Coefficient  Inverse  Problem  (CIP)  for  a  Partial  Differential  Equation  (PDE)  is  a  prob¬ 
lem  of  the  reconstruction  of  an  unknown  coefficient  of  that  PDE  from  the  boundary  measure¬ 
ments.  Because  of  many  dangers  on  the  battleheld,  the  Army  is  interested  in  the  minimal 
number  of  measurements.  Thus,  only  CIPs  with  a  single  measurement  event  are  considered 
in  this  project.  In  other  words,  either  a  single  location  of  the  point  source  or  a  single  direc¬ 
tion  of  the  initializing  plane  wave  is  considered.  It  was  observed  that  the  best  is  a  two-stage 
numerical  procedure: 

Stage  1.  The  globally  convergent  method  provides  a  good  approximation  for  the  solution 
without  any  a  priori  knowledge  of  a  small  neighborhood  of  this  solution. 

Stage  2.  The  locally  convergent  Adaptive  Finite  Element  Method  (adaptivity) 

The  hrst  stage  gives  us  accurate  values  of  dielectric  constants.  The  second  stage  provides 
accurate  images  of  shapes  of  targets.  It  is  important  that  the  adaptivity  uses  the  solution  of 
the  hrst  stage  as  its  starting  point. 


2  The  List  of  Main  Results  Of  This  Project 

The  hrst  paper  about  the  globally  convergent  method  of  this  paper  was  published  in  2008 
[25].  Since  then  the  authors  of  [25]  have  developed  this  method  much  further,  see  items  ^1-Q 
in  Introduction.  Our  group  has  pioneered  all  main  results  listed  below.  The  main  results  of 
this  project  are: 

1.  The  analytical  proof  of  the  global  convergence  property  using  a  sophisticated  mathe¬ 
matical  apparatus  [1,  16,  5]. 

2.  The  development  of  a  sophisticated  analytical  apparatus  for  establishing  the  relaxation 
property  of  the  adaptivity  technique  [1,  7]. 
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3.  Numerical  implementations  of  resulting  algorithms. 

4.  Numerical  verifications  of  resulting  algorithms  on  computationally  simulated  data  [1, 
16], 

5.  Assembling  an  experimental  apparatus  in  Microwave  Laboratory  of  University  of  North 
Carolina  at  Charlotte. 

6.  Verihcation  of  the  globally  convergent  numerical  method  on  backscattering  experimen¬ 
tal  data  for  targets  standing  in  air  [2,  3,  4,  6].  Targets  mimic  explosives. 

7.  Verihcation  of  the  globally  convergent  numerical  method  on  backscattering  experimen¬ 
tal  data  for  targets  buried  in  the  ground  [5].  This  case  is  much  more  complicated  than 
the  case  of  targets  in  air. 

8.  An  experimental  and  numerical  reconstruction  evidence  of  the  super  resolution  phe¬ 
nomenon  [5]. 

9.  Addressing  a  need  of  the  Army  via  successful  work  with  experimental  data  collected 
by  the  Forward  Looking  Radar  of  US  Army  Research  Laboratory  (ARL)  [17,  18,  19]. 
The  globally  convergent  method  of  this  project  was  used. 

10.  Transfer  of  a  ready-to-use  software  to  ARL  resulting  from  item  ^^9.  This  software 
works  with  the  real  data  of  the  Forward  Looking  Radar  of  ARL  [59]. 

11.  The  use  of  experimental  data  of  item  ^^9  for  a  comparison  of  performances  of  the 
globally  convergent  numerical  method  of  this  project  and  the  classical  Krein  eqnation 
method  [19].  It  was  established  that  while  the  hrst  method  works  well,  the  second  one 
fails  for  these  data. 

12.  Four  (4)  presentations  to  Mr.  Brian  Bnrns,  Drs.  Anders  Snllivan  and  Lam  Nguyen, 
ARL  engineers,  and  to  Dr.  Joseph  D.  Myers,  the  Program  Manager  of  the  Numerical 
Analysis  Program  of  ARO. 

13.  Three  joint  publications  with  ARL  engineers  Drs.  Anders  Snllivan  and  Lam  Nguyen 
[17,  18,  19]. 

14.  A  new  globally  convergent  nnmerical  method  based  on  the  Carleman  Weight  Fnnction 

[11]. 

15.  The  hrst  solution  of  a  long  standing  problem  about  uniqueness  of  a  phaseless  3-d 
inverse  problem  of  qnantnm  scattering  [8,  9,  10].  This  was  an  open  qnestion  since  the 
pnblication  of  the  well  known  book  [38]  in  1977. 

16.  Pnblication  of  a  book  summarizing  some  results  of  the  project  [1].  This  book  was 
published  by  Springer,  the  World  most  prestigious  publisher  of  the  scientihc  literature. 

17.  Publication  of  two  survey  papers  [7,  13]. 
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3  The  Theory  of  The  Globally  Convergent  Numerical 
Method 

3.1  Statement  of  the  CIP 

As  the  forward  problem,  we  consider  the  following  Cauchy  problem 

er{x)utt  =  Am  in  x  (0,  cxd)  ,  (1) 

u{x,0)  =  0,Ut{x,0)  =  6  {x  —  Xq)  .  (2) 

Here  er{x)  is  the  spatially  distributed  variable  dielectric  constant.  It  is  well  known  that 
equation  (1)  can  be  derived  from  Maxwell’s  equations  only  in  the  2-d  case.  Also,  even 
though  it  cannot  be  derived  in  the  3-d  case  for  er{x)  ^  const.,  it  was  successfully  used 
to  model  electric  waves  propagation  in  our  works  with  experimental  data  in  [1,  2,  3,  4,  5]. 
The  reason  of  the  success  with  this  equation,  instead  of  the  Maxwell’s  system,  for  real  data 
was  explained  numerically  in  [29].  Namely,  it  was  shown  in  [29]  that  if  the  incident  electric 
held  has  the  form  Einc  =  (0,i?2,0),  then  the  component  E2{x,t)  dominates  components 
El  {x,  t) ,  E^  {x,  t)  and  the  propagation  of  E2  [x,  t)  is  well  governed  by  equation  (1). 

Let  n  C  R^  be  a  convex  bounded  domain  with  the  boundary  dQ  G  C^.  We  assume  that 
the  coefficient  Er  {x)  of  equation  (1)  is  such  that 

Er  (x)  G  [1,  d]  ,  Er  {x)  =  1  for  X  G  R^\r2,  (3) 

Er{x)  G  ^^(R^),  (4) 

where  d  =  const.  >  1  characterized  the  upper  bound  of  the  function  Er  {x) .  Since  we  do  not 
impose  any  smallness  conditions  on  the  number  d  =  d  —  1,  then  our  results  are  not  “local” 
ones. 

Coefficient  Inverse  Problem  3.1  (CIP3.1).  Suppose  that  the  coefficient  Er  {x)  satis¬ 
fies  (3),  (4).  Assume  that  the  function  Er  {x)  is  unknown  in  the  domain  hi.  Determine  the 
function  Er  {x)  for  x  E  fl,  assuming  that  the  following  function  g  {x,  t)  is  known  for  a  single 
source  position  Xq  ^  kl 

u  {x,  t)  |aox(o,oo)=  g  {x,  t) .  (5) 

The  assumption  Er  (x)  =  1  for  x  G  R^\r2  means  that  one  has  air  outside  of  the  medium 
of  interest  hi.  The  inequality  Er  (x)  >  1  is  because  the  speed  of  EM  waves  propagation 

in  the  medium  is  less  than  one  in  the  air.  The  function  g  (x,  t)  models  time  dependent 

measurements  of  the  wave  field  at  the  boundary  of  the  domain  of  interest.  The  assumption 
that  the  function  g  [x,  t)  is  known  on  the  infinite  time  interval  t  G  (0,  cxd)  rather  than  on  a 
finite  one,  is  not  a  serious  restriction  from  the  computational  point  of  view.  Indeed,  we  work 
with  the  Laplace  transform  (6)  and  the  kernel  of  the  integral  (6)  decays  very  rapidly. 
Therefore,  in  fact  only  values  of  the  function  g{x,t)  for  t  G  (0,a)  with  a  small  number 
a  are  used  in  our  method.  Besides,  our  data  pre-processing  procedure  for  experimental 
data  [1,  2,  3,  4,  5]  shows  that  we  can  work  with  only  a  small  piece  of  the  time  dependent 
experimental  curves. 
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3.2  Integral  differential  equation 

Below  (7^+“  denote  Holder  spaces,  where  /c  >  0  is  an  integer  and  a  G  (0, 1) .  Consider  the 
Laplace  transform  of  the  functions  m, 


w{x,s)  =  J  u{x,t)e  for  s  >  s  =  const.  >  0,  (6) 

0 

where  s  is  a  sufficiently  large  number.  In  our  numerical  studies  we  choose  s  numerically. 
Then  (1),  (2)  and  Theorem  2.7.1  of  [1]  imply  that 

Aw  —  s'^Er  {x)  W  =  —6{x  —  Xo),XEM.^,ys>S,  (7) 

lim  w{x,s)  =  0,  \/s  >  s.  (8) 

|x|^oo 

Theorem  3.1  is  a  reformulation  of  Theorem  2.7.2  of  the  book  [1].  A  similar  theorem  for  the 
case  when  the  point  source  in  (2)  is  replaced  with  an  incident  plane  wave  was  proven  in  [5]. 

Theorem  3.1.  Let  the  souree  Xq  and  the  funetion  Sr  {x)  satisfies  conditions  (3)  and 
also  Er  G  C"  (M^) .  Let  wi  {x,  s)  and  Wd  {x,  s)  be  solutions  of  the  problem  (7),  (8)  for  =  1 
and  Er  =  d  respectively, 

exp  (— s  |a;  —  Xol)  exp  (—sy/d\x  —  Xo\) 

[X:  S)  =  - — - ^ - ,  Wd  [x,  s)  =  - — - ^ - . 

Ati\X  —  Xq\  Ati\X  —  Xq\ 

Then  for  any  s  >  0  there  exists  unique  solution  w{x,s)  of  the  problem  (7),  (8),  which  is 
represented  in  the  form 

w  {x,  s)  =  wi  {x,  s)  +w  {x,s)  ,w  e  (M^)  .  (9) 

Furthermore, 

Wd  {x,  s)  <  W  {x,  s)  <  Wi  {x,  s)  ,  Vx  7^  Xq.  (10) 

Also,  the  Laplace  transform  (6)  w{x,s)  of  the  function  u{x,t)  satisfies  (9),  (10)  for  suffi¬ 
ciently  large  .s  and  s  >  s. 

Thus,  below  we  consider  only  those  solutions  of  the  problem  (7),  (8),  which  satisfy  con¬ 
ditions  (8),  (10).  Since  by  (10)  tc  >  0,  then  we  can  consider  the  function  v  =  Intc/s^.  Hence, 
recalling  that  Xq  ^  H,  we  obtain  that  (7)  leads  to 

Av  -f  s^  |Vn|^  =  Er{x),  X  E  Ll,  (11) 

v\dn  =  9?(a;,s) ,  s  G  [s,s] ,  (12) 

where  the  function  <p{x,s)  is  generated  by  the  function  g{x,t)  in  (5).  Introduce  a  new 
function  q  [x,  s)  =  dgV  {x,  s) .  Then  under  certain  non-restrictive  conditions 


D“(v)  =  O  I  i)  ,  CJ (g)  =  O  (1  I  ,  s  -»  oo;  |a|  <  2, 


(13) 


=  —  q  {x,t)  dr. 


V  {x,  s) 


(14) 


We  represent  the  integral  (14)  as 


V  {x,  s)  =  —  /  q  {x,  t)  dr  +  V  {x,  s) , 


(15) 


where  s  >  s  is  a  large  parameter  which  should  be  chosen  in  numerical  experiments.  Actually, 
s  is  one  of  regularization  parameters  of  our  method.  We  call  V  {x,  s)  the  tail  function, 


V  {x,  s)  =  —  /  q  {x,  t)  dr. 


By  (15) 
By  (13) 


Tw  lnt(;(a;,  s)  ^ 

V  (x,  s)  =  — =2 — 1  x  e 


=  O 


-fc+i 


,  /c  =  0, 1;  s  — )■  cxD. 


We  obtain  from  (11)-(17)  the  following  nonlinear  integral  differential  equation 

r  \  r  l' 

,2, 


s 

\  1  1 

/  Vg  {x,  r)  dr  +  2s 

<i 

Si¬ 

's 

kJ 

s 

L.  J 

+  2s^VqW  —  2s'VV  ■  j  Vg  {x,  r)  dr  +  2s  (VB)^  =  0. 
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Let  Ip  {x,  s)  =  dgip  {x,  s) .  Then  (12)  implies  that 


(16) 

(17) 


(18) 


Slasi  =  V’ (i:,  s) ,  s  6  |s,  s]  ■  (19) 

Hence,  we  need  to  solve  the  problem  (18),  (19).  Here  the  truncation  parameter  s  is 
the  regularization  parameter  of  our  numerical  method.  The  presence  of  integrals  in  (18) 
implies  the  nonlinearity,  which  is  the  main  difficulty  here.  If  both  functions  q  and  V  are 
approximated  well  together  with  their  x— derivatives  up  to  the  second  order,  then  the  target 
unknown  coefficient  Er  {x)  can  also  be  approximated  well  via  (11),  where  the  function  v  is 
computed  via  (15).  Equation  (18)  contains  two  unknown  functions  q  and  V.  The  reason  why 
we  can  approximate  both  of  them  is  that  we  treat  them  differently:  while  we  approximate 
the  function  q  via  outer  iterations,  the  function  V  is  approximated  via  inner  iterations.  Thus, 
below  we  focus  on  the  following  question:  How  to  solve  the  problem  (18),  (19)  numerically? 
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3.3  Layer  stripping  with  respect  to  s 

We  approximate  the  function  q  [x,  s)  as  a  piecewise  constant  function  with  respect  to  the 
pseudo  frequency  s.  That  is,  we  assume  consider  a  partition  of  the  interval  [s,  s]  with  the 
sufficiently  small  grid  step  size  h, 

S_  =  Sjv  ^  ^N—l  Si  Sq  =  S ^  h  =  Si— I  —  Si- 

We  assume  that  q  {x,  s)  =  qn  {x)  for  s  G  (s„,  .  The  boundary  condition  (19)  is  approxi¬ 

mated  as 

qnlan  =  (20) 

where  is  the  average  of  the  function  ip  over  the  interval  (s„,  s„_i) .  Rewrite  (18)  for 
s  G  {sn,  s„_i]  using  this  piecewise  constant  approximation.  Then  multiply  the  resulting 
approximate  equation  by  the  s-dependent  Carleman  Weight  Function  (CWF)  of  the  form 

^n,ii  ("S)  exp  [  jj,  |s  ,  s  G  (s^,  ,  (21) 

and  integrate  with  respect  to  s  G  (s„,s„_i] .  Here  /i  >>  1  is  a  large  parameter  of  ones 
choice.  Usually  we  choose  /i  =  50.  We  obtain  the  following  approximate  equation  in  H  for 
the  function  (x)  ,n  =  1, ...,  N 


n—1 


Ln  {Qti)  •  -^In  j  ^  ^  Qj  ^  Qn 


(22) 


'’n—l 


n—1 


Bn  (V^n)  ~  ^2,nh^  |  -|-  2H2,nVUn  |  h  —  H2,n  (VRi 


d=0 


j=0 


We  have  intentionally  inserted  dependence  of  the  tail  function  14  from  the  iteration  number 
n  here  because  we  will  approximate  these  functions  iteratively.  In  (22)  Ai^n  =  Ri,n  (yW,  h) , 
A2,n  =  R2,n  (p-,  h)  and  =  Bn  {fi,  h)  are  certain  numbers  depending  on  jj,  and  h,  see  specihc 
formulas  in  (22).  It  is  convenient  to  set  in  (22) 


go  =  0. 

Since  boundary  value  problems  (20),  (22)  are  actually  generated  by  equation  (18),  which 
contains  Volterra-like  s-integrals,  then  these  problems  can  be  solved  sequentially  starting 
from  qi.  As  to  (22),  an  important  point  is  that  \Bn  (/i,  h)\  <  Ss"^/ fx  for  jxh  >  1  (22).  We  have 
used  yU  =  50  in  our  computations.  Hence,  assuming  that  fx  »  1,  we  ignore  the  nonlinear 
term  in  (22)  via  setting 

Bn  (Vgn)'  :=  0.  (23) 

This  allows  us  to  solve  a  linear  problem  for  each 
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3.4  The  algorithm 


On  each  iterative  step  n  we  approximate  both  the  function  and  the  tail  function  14,  see 
Remark  5.1.  Each  iterative  step  requires  an  approximate  solution  of  the  boundary  value 
problem  (20),  (22).  First,  we  choose  an  initial  tail  function  G  (O).  The  hrst 

choice  Vi  1  for  the  tail  function  is  described  in  subsection  3.6.  For  each  we  have  inner 
iterations  to  update  tails.  On  these  iterations  we  compute  functions  g„^fc,  /c  =  1,  ...,m.  Here 
the  number  of  inner  iterations  m  is  chosen  in  numerical  experiments.  The  criterion  for  the 
choice  is  the  stabilization  of  iteratively  computed  tail  functions  Vn,k  at  the  iteration  k  \=  m. 

Step  n^,  where  n, /c  >  1.  Recall  that  go  =  0.  Suppose  that  functions  qj  G  (O)  ,j  = 
l,...,n  —  1  and  tails  Ri, ...,  14_i,  14,^  G  (O)  are  constructed.  To  construct  the  function 
<in,ki  we  use  the  FEM  to  solve  the  following  Dirichlet  boundary  value  problem  for  the  elliptic 
equation  in  H 


Ag, 


n,/c 


V14,fc  Vq, 


n,k  — 


-^2,nh^ 


n—1 


+  2H2^nV14,fc  ■  j  h  Vgj  —  {S^n,k) 


j=0 


(24) 


qn,k\dn  = 

To  reconstruct  an  approximation  (x)  for  the  function  Sr  (x) ,  we  hrst  use  the  following 
discrete  analog  of  (15) 


n—1 

Vn,k  {x,  Sn)  =  -hqn,k  (x)  -  qj  (x)  +  14, fc  (x) . 

1=1 


Next,  let  C  H  be  a  subdomain  of  such  that 


(25) 


Q'  c  H,  dQ'  ndQ  =  0. 


(26) 


Consider  a  function  y  (x) , 


f  1  in  H', 

X  G  (M^)  ,x(x)  =  between  0  and  1  in  Q\Q', 

0  in 

(27) 

For  any  function  a 

G  C"  (fl)  dehne 

a  (x)  =  (1  —  X  (x))  +  X  (x)  a  (x) ,  X  G  M^. 

(28) 

Hence, 

f  l,x  E  R^\H, 

a  (x)  =  <  a  (x)  in  H', 

y  >  1  in  H,  if  a  (x)  >  1  in  H. 

(29) 

11 


Next,  using  (11),  (25),  (28)  and  (29),  we  set 


(30) 

(31) 


fn,k{x)  :  =  AVn,k{x,Sn)  +  sl\VVn,k{x,Sn)f  ,n>l,X  en, 

(x)  :  =  fn,k  (a:)  • 


Next,  we  solve  the  forward  problem  (7),  (8)  with  {x)  :=  {x)  ,s:=s  and  obtain  the 

function  Wn,k  {x,s)  satisfying  conditions  (9),  (10)  (Theorem  .1).  Next,  using  (16),  we  set  for 
the  new  tail 


hn,fc+i  (a^) 


In  Wr 


,k  [X,s) 


We  continue  these  iterations  with  respect  to  tails  for  k 


,  a;  G 

=  1, ...,  m.Next,  we  set 


,(n) 


x)  := 


x') ,  Qyi  (a^)  •  Qn,m  (a^)  ,  Vji  (x)  .  h7i,m  (a^) 


14i+i,i  (a^)  for  X  E  Q. 


We  stop  iterations  with  respect  to  n  at  n  :=  N  E  [1,1V] ,  where  iV  is  a  certain  number 
at  which  convergence  occurs.  Stopping  criteria  are  discussed  in  detail  in  our  publications 
[1,  16,  2,  3,  4,  5]and  they  go  along  well  with  one  of  backbone  ideas  of  the  theory  of  ill-posed 
problems.  By  this  idea,  the  iteration  number  can  be  chosen  as  a  regularization  parameter, 
see  pages  156  and  157  of  the  book  [40].  Furthermore,  we  have  constantly  observed  that 
results  of  our  numerical  studies  go  along  well  with  convergence  criteria. 


3.5  Estimates  of  tail  functions 

The  most  difficult  question  on  the  way  of  proving  our  advanced  convergence  theorem  was  the 
question  about  estimating  gradients  of  tail  functions  in  Holder  norms  This  was  done 

in  Theorem  2. 9. 1.2  of  the  book  [1]  as  well  as  in  Theorem  4.2  of  [16].  We  now  reformulate 
these  theorems. 

Let  Hi  C  be  a  hnite  domain  with  dkli  E  and  such  that  H  C  Hi,  (9Hn(9Hi  =  0.  Let 
d  >  1  be  the  number  from  (3)  and  d  >  d  be  another  number.  Introduce  the  set  of  functions 
P  (d,  d)  as 

P  (d,  d)  =  jsr  {x)  E  C“  (H)  :  <  d  +  1,1  <  Er  {x)  <  d  -|-  1,  Vx  G  h|  . 

For  each  function  Er  E  P  (d,  d)  we  construct  the  function  Er  (x)  by  the  formula  (28),  where 
the  function  x{x)  is  dehned  in  (27).  Let  {x,s)  be  the  solution  of  the  problem  (7)-(10) 
with  Er  :=  '£r,s  :=  s.  Using  (16),  dehne 


Vsr  (X,  S)  = 


In 


wg 


(32) 


Following  again  one  of  the  main  concepts  of  the  theory  of  Ill- Posed  Problems  [1] ,  we  assume 
in  Theorem  3.2  the  existence  of  the  exact  solution  e*  (x)  for  the  unperturbed  exact  data 
g*  (x,t)  in  (3). 
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Theorem  3.2.  Let  the  parameter  s  >  1  and  xq  ^  f2i.  Let  e*  (x)  G  P[d,d)  be  the 
exaet  solution  of  CIP3.1  with  the  noiseless  data  g*  {x,t)  in  (3)  and  V^*  {x,s)  be  the  exaet 
tail  defined  by  the  formula  (32).  Then  there  exists  a  constant  B  =  B  (f2,  f2i,  x,  a;o,  s,  Xq)  >  2 
depending  only  on  listed  parameters  such  that  the  following  estimates  hold 

II  <  -B,  <  B  ,\/€r  E  P  (^d,  d^  , 

II ^  E  P  (ci,  d)j . 


3.6  The  first  tail 


In  the  description  of  the  algorithm  in  subsection  3.4  we  have  left  open  the  question  about 
the  choice  of  the  hrst  tail  function  lh,i-  This  question  is  addressed  in  this  subsection.  By 
(16)  the  asymptotic  behavior  of  the  exact  tail  is 


^  +  O  ,s  ^  oo,x  E  il. 

s  V-s  / 


for  a  certain  function  p*  [x) .  We  truncate  the  second  term  of  this  asymptotic  behavior.  Thus, 
our  Approximate  Mathematical  Model  consists  of  the  following  assumption. 

Assumption.  There  exists  a  function  p*  (x)  E  (hi)  such  that  the  exact  tail  function 
V*  {x,  s)  has  the  form 

Ve*{x,s):='^^-^,ys>s.  (33) 

Furthermore,  by  (32)  we  assume  that 


p*  {x) 
s 


lnM;*(a;,s)  ^  _ 

- ^ - ,  Vs  >  s. 

qZ 


Since  q*  {x,  s)  =  dglfi*  {x,  s)  for  s  >  s,  we  derive  from  (33)  that 


(34) 


q‘(x,s)  =  -^.  (35) 

s 

Substituting  (33)  and  (35)  in  the  equation  (18)  and  boundary  data  (19)  for  the  function 
q*  {x,  s)  an  setting  there  s  :=  s,  we  obtain  the  following  approximate  Dirichlet  boundary 
value  problem  for  the  function  p*  [x) 


Ap*  =  OinQ,  p*  E  (^2+“  (fi)  ,  (36) 

P*\dn  =  -s^fj*ix,s),  (37) 


where  fj*  (x,  s)  is  the  exact  function  tjj  (x,  s) ,  which  corresponds  to  the  function  g*  {x,  t) .  The 
approximate  equation  (36)  is  valid  only  within  the  framework  of  Assumption.  Although  this 
equation  is  linear,  formula  (11)  for  the  reconstruction  of  the  target  coefficient  is  nonlinear. 
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Recall  that  by  (19)  q  {x,  s)  =  ip  (x,  s) ,  V  {x,  s)  G  (911  x  [s,  s]  .  Assume  that 

tp  (X,  S)  e  C'2+a  ^  g  _  (-33^ 

Consider  the  solution  p  (x)  of  the  following  boundary  value  problem 

Ap  =  0  in  n,  p  G  (n)  ,  (39) 

plan  =  -s^P’(x,s) .  (40) 

Given  (38),  the  Schauder  theorem  implies  that  there  exists  unique  solution  p  of  the  problem 
(39),  (40). 

As  the  hrst  guess  for  the  tail  function  we  take 


n,.  (X)  :=  dh. 

s 

It  follows  from  (33)- (41)  that 

l|Ci,i  -  <  M  \\ip  {x,s)  -  Ip*  (x,s)||^2+.(af^) 


(41) 


(42) 


where  M  =  M  (H)  =  const.  >  0. 

Remarks: 

1.  We  point  out  that  the  truncation  of  the  asymptotic  series  with  respect  to  1/s  is  done 
in  Assumption  only  on  the  first  iteration,  i.e.  only  for  the  hrst  tail.  No  assumptions  are 
imposed  on  follow  up  iterations. 

2.  An  interesting  conclusion  which  follows  immediately  from  (42)  is  that  we  obtain  a 
good  approximation  for  the  solution  of  our  CIP  already  from  the  hrst  tail.  In  other  words, 
a  good  approximation  is  obtained  already  on  the  first  iteration  of  our  method.  And  we 
have  constantly  observed  this  computationally.  The  error  of  this  approximation  depends 
only  on  the  error  in  the  boundary  data.  It  follows  from  here  and  from  (11),  (15)  that  we 
obtain  a  good  approximation  for  the  target  coefficient  already  from  the  first  tail.  Theorem 
3.3  guarantees  that  all  other  solutions  obtained  in  the  above  iterative  process  also  provide 
good  approximations,  as  long  as  the  number  of  iterations  is  not  too  large.  This  means 
that  we  should  develop  numerically  a  stopping  criterion  to  stop  iterations.  Suppose  now 
that  iterations  are  stopped  before  this  stopping  criterion  is  met.  In  particular,  they  can  be 
stopped  just  on  the  hrst  iteration.  In  this  case  we  can  apply  the  second  stage  of  our  two-stage 
numerical  procedure  [1].  Namely,  we  could  apply  a  locally  convergent  adaptivity  technique 
to  rehne  solution.  We  would  take  the  solution  obtained  on  the  globally  convergent  stage  as 
the  starting  point  of  iterations. 


3.7  Global  convergence  theorem 

The  proof  of  Theorem  3.3  was  hrst  published  in  our  book  [1]:  see  Theorem  2.9.4  there.  Next, 
it  was  published  in  [16].  Theorem  3.3  was  extended  in  [5]  to  the  case  when  the  point  source 


14 


is  replaced  by  an  incident  plane  wave.  The  proof  of  Theorem  3.3  essentially  uses  Theorem 
3.2  and  Assumption. 

Let  G  dfl  be  the  boundary  data  (20)  for  the  functions  qn{x)  generated  by 

the  data  function  'ip{x),x  G  dfl  in  (19).  Let  'ip*{x),x  G  dfl  be  the  function  ■^(x)  which 
corresponds  to  the  exact  solution  e*{x),xEfl  of  CIP3.1,  which  was  introduced  in  Theorem 
3.2.  Let  'ipn{^)jX  G  <911  be  corresponding  exact  functions  Let  the  small  number 

a  G  (0, 1)  represents  the  level  of  the  error  in  the  data  iIj{x).  Then  it  is  natural  to  assume 
that 


c‘^+°‘(an) 


<  C*  {h  +  a) ,  C*  =  const.  >  1. 


(43) 


Theorem  3.3.  Let  Assumption  be  valid.  Consider  the  algorithm  of  subsection  3.f.  Let 
the  parameter  s  >  1  and  xq  ^  fli.  Let  e*  (x)  G  P  (d,  d)  be  the  exact  solution  of  CIP3.1  with 
the  noiseless  data  g*  {x,t)  in  (3)  and  Vg*  (x,s)  be  the  exact  tail  defined  by  the  formula  (32). 
Assume  that  (43)  holds.  Introduce  the  error  parameter  rj, 


T]  =  2{h  +  a) . 


Assume  that  the  parameter  p,  in  the  Carleman  Weight  Function  (21)  is  so  large  that 


L  > 


8  (sC*)^ 
V 


Also,  assume  that  in  that  algorithm  all  functions  {x)  >  l,a;  G  fl.  Then  there  exists  a 
constant  B  =  B  (hi,  hli,  y,  xq,  s,  Xq)  >  2  depending  only  on  listed  parameters  such  that  if  the 
parameter  rj  is  so  small  that 


T]  G  (0,r/o)  where 


(44) 


then  the  following  estimate  holds  with  a  number  u  =  u  {B,  N,  m)  G  (0, 1)  depending  only  on 
listed  parameters 

-<|lc«(o)  ^  ^  (O’^o)  • 

In  other  words,  the  algorithm  of  subsection  3.4  is  globally  convergent. 

In  particular,  (44)  requires  that  the  error  parameter  g  and  the  total  number  of  iterations 
Nm  should  be  connected  with  each  other.  We  again  refer  to  pages  156  and  157  of  the 
classical  book  about  the  theory  of  Ill- Posed  Problems  [40].  It  follows  from  there  that  such 
connection  should  be  in  place  even  for  simpler  ill-posed  problems. 


4  The  Theory  Of  the  Image  Refinement  via  Adaptive 
Finite  Element  Method  (Adaptivity) 

As  it  was  pointed  out  in  Introduction,  we  have  constantly  observed  that  the  best  approach 
is  to  apply  a  two-stage  numerical  procedure.  On  the  hrst  stage  the  globally  convergent 
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numerical  method  delivers  a  function  which  is  sufficiently  close  to  the  exact  coefficient.  In 
fact,  the  first  stage  provides  accurate  locations  and  accurate  values  of  refractive  indices 
n  {x)  =  of  targets  mimicking  explosives.  On  the  second  stage  the  locally  convergent 

method,  the  adaptivity,  is  applied  to  rehne  images.  The  adaptivity  uses  the  solution  of  the 
hrst  stage  as  its  starting  point.  The  adaptivity  provides  accurate  shapes  of  those  targets. 
Therefore,  the  two-stage  procedure  provides  all  three  components  of  interest  of  targets: 

1.  Refractive  indices. 

2.  Locations. 

3.  Shapes. 

The  classical  Tikhonov  functional  for  solving  ill-posed  problems,  including  CIPs  is  known 
for  a  long  time,  see,  e.g.  the  books  [1,  40].  However,  what  makes  the  adaptivity  attractive 
is  that  it  minimizes  this  functional  on  a  sequence  of  locally  rehned  meshes  of  hnite  elements 
rather  than  on  a  single  mesh.  Although  the  adaptivity  for  ill-posed  problems  has  been  known 
since  2001,  see,  e.g.  [30]- [35],  a  rigorous  proof  of  the  key  fact  that  mesh  rehnements  indeed 
provide  a  better  accuracy  was  absent.  We  call  the  latter  property  relaxation.  The  major 
obstacle  for  the  proof  of  the  relaxation  property  is  the  ill-posed  nature  of  inverse  problems. 
The  first  rigorous  proof  of  this  fact  was  obtained  in  2010  [28].  Next,  this  result  was  rehned 
in  the  book  [1].  A  survey  was  published  in  [7]. 

In  this  section  we  present  our  results  about  the  relaxation  of  the  adaptivity.  We  point 
out  that  it  was  natural  to  prove  hrst  the  relaxation  property  for  an  abstract  operator  F. 
However,  it  is  well  known  in  the  theory  of  ill-posed  problems  that  it  is  usually  not  easy 
to  transform  this  result  to  the  case  of  our  specihc  CIP  3.1.  Nevertheless,  this  was  done  in 
[1,  7,  28]. 

4.1  The  space  of  finite  elements 

As  the  hrst  step  for  establishing  the  relaxation  property  of  the  adaptivity,  we  introduce  the 
space  of  hnite  elements.  It  is  worthy  to  point  out  that  the  group  led  by  the  PI  was  the  first 
one  who  has  introduced  this  space  in  2010  [28]. 

Let  H  C  M"',n  =  2,3  be  a  bounded  domain.  Consider  a  discretization  of  hi  by  an 
unstructured  mesh  T  using  non-overlapping  elements  K.  For  the  elements  K  are  triangles 
or  quadrilaterals  and  for  R^  tetrahedrons  or  hexahedrons  such  that  T  =  Ki, Ki,  where  I 
is  the  number  of  elements  in  hi,  and 

D  =  UkgtK  =  Ki  U  K2...  U  Ki. 

We  obtain  a  polygonal  domain  D  and  assume  for  brevity  that  D  =  Q. 

Following  section  76.4  of  the  book  [42],  consider  piecewise  linear  functions  {ej  {x,  C 

C  (H),  which  are  called  test  functions.  Functions  {ej  (a;,  T)}^^  are  linearly  independent  in 
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f2.  Here,  N  is  the  total  number  of  nodes  in  the  mesh  T.  Let  {iVj}  be  the  set  of  nodal  points 
of  triangle/tetrahedra  K  for  all  K  eT.  Then 

We  introduce  the  hnite  element  space  Vh  as 

Vh  =  {n(a;)  G  H\n)  :  v  E  0(11),  v\k  G  Pi{K)  ViL  G  T}, 

where  Pi{K)  denotes  the  set  of  piecewise-linear  functions  on  K.  The  hnite  dimensional  hnite 
element  space  14  is  constructed  such  that  14  C  V. 

Let  hx  be  the  diameter  of  element  K  which  we  dehne  as  the  longest  side  of  K  and  r 
be  the  radius  of  the  maximal  circle/sphere  inscribed  in  K.  We  impose  the  shape  regularity 
assumption  for  all  triangles/tetrahedra  uniformly  for  all  possible  triangulations  T  which  we 
consider.  Specihcally,  we  assume  that 


4  hif  4  'ra2,  ai,  02  =  const.  >  0,  \/K  E  T,  T,  (45) 

where  numbers  Oi,  02  are  independent  on  the  triangulation  T.  Let  hmax  (T)  and  hmin  (T)  be 
respectively  the  maximal  and  minimal  diameters  of  triangles/tetrahedra  of  the  triangulation 
T.  We  assume  everywhere  below  that 


hmin  (T) 
hmax  (T) 


<  ct,^T 


(46) 


for  a  certain  positive  constant  ct-  Obviously,  the  number  of  all  possible  triangulations 
satisfying  (45),  (46)  is  hnite.  Thus,  we  introduce  the  following  hnite  dimensional  linear 
space  H, 

H  =  \JVhiT),  VT  satisfying  (45),  (46). 

T 

Hence, 

dim  H  <  00,  H  <Z  {C  (O)  fl  P[^  (O))  ,  dxj  E  (O) ,  yfEH.  (47) 

In  (47)  ”c”  means  the  inclusion  of  sets.  We  equip  H  with  the  same  inner  product  as 
the  one  in  L2  (O) .  Denote  (, )  and  H-H  the  inner  product  and  the  norm  in  H  respectively, 
Wflln  II/IIl2(o)  ll/ll  Everywhere  in  section  4  77  is  this  space.  We  view  the 

space  H  as  an  “ideal”  space  of  very  hue  hnite  elements,  which  cannot  be  reached  in  practical 
computations.  At  the  same  time,  all  other  spaces  of  hnite  elements  we  work  with  below  are 
subspaces  of  H.  In  particular,  this  means  that  we  assume  without  further  mentioning  that 
(45)  and  (46)  are  valid  for  all  meshes  considered  below. 

Keeping  in  mind  the  mesh  rehnement  process  in  the  adaptivity,  we  now  explain  how  do 
we  construct  triangulations  {T„}  as  well  as  corresponding  subspaces  {M„}  of  the  space  H 
which  correspond  to  mesh  rehnements.  Consider  the  hrst  triangulation  Ti  with  rather  coarse 
mesh.  We  set  Mi  :=  14  (Ti)  C  H.  Suppose  that  the  pair  (Tn,Mn)  is  constructed  after  n 
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mesh  refinements  and  that  the  basis  functions  in  the  space  are  {ej  ■  We  now 

want  to  refine  the  mesh  again.  We  define  the  pair  (T^+i,  as  follows.  We  refine  the 

mesh  in  the  standard  manner  as  it  is  usually  done  when  working  with  triangular /tetrahedron 
finite  elements.  When  doing  so,  we  keep  (45).  Hence,  we  obtain  both  the  triangulation  T„+i 
and  the  corresponding  test  functions  {ej  (x,  T„+i)}^/h  It  is  well  known  that  test  functions 

{cj  are  linearly  dependent  from  new  test  functions  {ej  (a:,  .  Thus,  we 

define  the  subspace  Mn+i  as 

Mn+i  :=  Span  (^{ej  (x,  j  . 

Therefore,  we  have  obtained  a  finite  set  of  linear  subspaces  of  the  space  H.  Each 

subspace  M„  corresponds  to  the  mesh  refinement  number  n,  Mn+i\Mn  7^  0  and 

Mn  C  Mn+i  C  H,ne[l,N  -1]. 

Let  /  be  the  identity  operator  on  H .  For  any  subspace  M  C  H,  let  Pm  ;  iL  — )■  M  be 
the  orthogonal  projection  operator  of  the  space  H  onto  its  subspace  M.  Denote  for  brevity 
Pn  :=  Pm„-  Let  hn  be  the  maximal  grid  step  size  of  T^.  Hence,  hn+i  ^  hn-  Let  //  be  the 
standard  interpolant  of  the  function  f  E  H  on  triangles/tetrahedra  of  T„,  see  section  76.4  of 
[42].  It  can  be  easily  derived  from  formula  (76.3)  of  [42]  that 

11/  -  /„'||  <  K  II  A„,  V/  €  H,  (48) 

where  K  =  K  (f2,r,  01,02)  =  const.  >  0.  Since  //  E  H,\/f  E  H,  then  by  one  of  well  known 
properties  of  orthogonal  projection  operators, 

ll/--Pn/ll  <  ||/-/,(||.  v/eff.  (49) 

Hence,  (48)  and  (49)  imply  that  with  a  different  constant  K  =  K  (f2,r,  01,02)  >  0 

11/  -  P„/||  s:  K  l|V/|7^,„|  fc„,V/  €  H.  (50) 

Since  77  is  a  finite  dimensional  space  in  which  all  norms  are  equivalent,  it  is  convenient  for 
us  to  rewrite  (50)  with  a  different  constant  K  =  K  (D,  r, ,  oi,  02)  >  0  as 


||a;  —  Pnx\\  ^  K  ||a;||  hn,  Vx  E  H.  (51) 

For  any  o  >  0  and  for  any  x  E  H  denote  Va{x)  =  {z  E  P[  ||a;  —  z\\  <  0}  . 

4.2  Relaxation 

Let  H2  be  another  Hilbert  space.  Let  G  C  77  be  an  open  bounded  set  and  F  :  G  ^  H2  he 
an  operator  which  is  continues,  has  the  Frechet  derivative  F'  {x)  ,'ix  E  G  and  the  operator 
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F'  [x)  is  Lipschitz  continous  on  the  set  G.  For  convenience,  we  assume  below  that  F  is 
one-to-one.  Consider  the  equation 


F  (x)  =  y,x  e  G.  (52) 

As  it  is  usually  done  in  the  regularization  theory  [1,  40],  we  assume  that  the  right  hand  side 
of  equation  (52)  is  given  with  a  small  error  S  G  (0, 1).  We  also  assume  that  there  exists  an 
“ideal”  exact  solution  x*  of  (52)  with  the  “ideal”  exact  data  y*.  Thus,  we  assume  that 

F{x*)=y*,x*  EG,\\y-y*\\^<5.  (53) 

Let  Xq  E  H  he  a  hrst  guess  for  the  exact  solution  x*.  Consider  the  Tikhonov  functional 
Ja  (x) , 

(x)  =  -\\F  {x)  -  y\\l  +  ^\\x  -  Xo\\^  ,x  E  G,  Xo  E  G,  (54) 

where  a  E  (0, 1)  is  the  regularization  parameter.  We  impose  a  rather  conventional  assump¬ 
tion  that 

a  =  a  (5)  =  (5^^,  y  =  const.  E  ^0,  .  (55) 

A  minimizer  Xa{5)  of  this  functional  is  called  regularized  solution.  Since  by  (47)  dim  77  <  oo, 
then  the  following  lemma  follows  immediately  from  the  Weierstrass  theorem. 

Lemma  4.1.  Let  F  be  the  operator  defined  above  in  this  section.  Then  there  exists  a 
regularized  solution  Xa{5)  E  G, 

inf  Ja  {x)  =  nun  Jq  {x)  =  Ja  {xq.)  .  (56) 

G  G 

Theorem  4.1  specihes  the  location  of  the  minimizer  Xa{s)  of  the  functional  (54). 
Theorem  4.1.  Let  Vi  (x*)  C  G.  Let  in  (54)  the  first  guess  xq  for  the  exact  solution  x* 
be  so  accurate  that 

||a;o-a;*||  <— .  (57) 

Then  there  exists  a  sufficiently  small  number  5o  =  5o  (-^i,  -^2,  y)  G  (0, 1)  such  that  for  every 
6  E  (0,  (5o)  and  for  a  =  a  {6)  satisfying  (55)  there  exists  unigue  regularized  solution  Xa{s) 
of  eguation  (52)  on  the  set  G.  Furthermore,  Xa{5)  E  V^3^i|.^{x*) .  In  addition,  the  gradient 
method  of  the  minimization  of  the  functional  Ja{5)  {x) ,  which  starts  at  xo,  converges  to 
Xa{5).  In  the  noiseless  case  with  <5  =  0  one  should  replace  “(5o  =  So{F,y)  E  (0,1)”  with 
«o  =  tto  {F)  G  (0, 1)  to  be  sufficiently  small  and  also  require  that  a  E  (0,  oq)  • 

Since  we  sequentially  minimize  the  Tikhonov  functional  on  subspaces  in  the 

adaptivity  procedure,  then  we  need  to  establish  hrst  the  existence  of  a  minimizer  on  each 
of  these  subspaces.  Theorem  4.2  ensures  both  existence  and  uniqueness  of  the  minimizer  of 
the  functional  Ja  on  each  subspace  of  the  space  77,  as  long  as  the  maximal  grid  step  size  of 
hnite  elements,  which  are  involved  in  that  subspace,  is  sufficiently  small. 
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Theorem  4.2.  Let  conditions  of  Theorem  4-1  hold.  Let  M  C  H  be  a  subspace  of  H. 
Assume  that  ||a;*||  <  -B,  where  the  number  B  >  0  is  known  in  advance.  Suppose  that  the 
maximal  grid  step  size  h  of  finite  elements  of  M  is  so  small  that 


h  < 


5BRK' 


(58) 


where  K  is  the  constant  in  (51)  and  the  constant  R  >  0  depends  only  on  the  operator 
F.  Furthermore,  assume  that  the  first  guess  xq  for  the  exact  solution  x*  in  the  functional 
Ja{s)  is  SO  accurate  that  (51)  is  in  place.  Then  there  exists  a  sufficiently  small  number 
<^o  =  <^o  {F,  fi)  G  (0, 1)  such  that  for  every  5  G  (0,(5o)  there  exists  unigue  minimizer  XM,a{S)  ^ 
G  0  M  of  the  functional  Jq,  on  the  set  GAM.  Furthermore,  XM,a{5)  ^  {x*)  fl  M.  Let 

Xa(5)  £  [x*)  be  the  regularized  solution  of  eguation  (52),  which  is  guaranteed  by  Theorem 

4.1.  Then  the  following  a  posteriori  error  estimate  holds 


We  are  now  ready  to  formulate  our  main  result  about  the  adaptivity,  which  is  the  relax¬ 
ation  Theorem  4.3 

Theorem  4.3  (relaxation).  Let  Mn  C  H  be  the  subspace  obtained  after  n  mesh  refine¬ 
ments,  as  described  in  subsection  4-L  Let  hn  be  the  maximal  grid  step  size  of  the  subspace 
M„.  Suppose  that  all  conditions  of  Theorem  4-2  hold  with  the  only  exception  that  the  subspace 
M  is  replaced  with  and  the  ineguality  (58)  is  replaced  with 


h„  < - . 

-  5BRK 

Let  S  G  (0,  (5o) ,  where  the  number  (5o  G  (0, 1)  is  defined  in  Theorem  4-2.  Let  Xn  G  {x*)  fl 
Mn  be  the  unigue  minimizer  of  the  functional  Ja  {x)  in  (54)  on  the  set  G  fl  (Theorem 
4-2).  Let  Xa{5)  G  Vg3ii/^{x*)  be  the  unigue  regularized  solution  (Theorem  4-k).  Assume  that 


7^  ^a{S)  •) 

i.  e.  Xa{s)  ^  Mn  ,  meaning  that  the  regularized  solution  is  not  yet  reached  after  n  mesh  refine¬ 
ments.  Let  7]  G  (0, 1).  Then  one  can  choose  the  maximal  grid  size  hn+i  =  hn+i  {F,  S,  B,  K,  rj)  G 
(0,  hn]  of  the  mesh  refinement  number  {n  +  1)  so  small  that 

ll^^n+l  Xa(S)  II  ^  ^  ll^^n  Xa(S)  ||  ) 

where  Xn+i  G  (x*)  fl  M^+i  is  the  unigue  minimizer  of  the  functional  (54)  on  the  set 
G  n  Mn+i .  Hence, 

ll^^n+l  Xa(5)  II  ^  ^  11^1  ^a(5)  ||  • 

Theorems  4. 1-4.3  are  formulated  for  an  abstract  operator  F.  They  need  to  be  specihed 
for  the  case  of  CIP3.1  (section  3.1),  which  is  a  non-trivial  task.  This  specihcation  was  done 
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in  [1,  7,  28].  Let  Ea{5)  (^r)  be  the  analog  of  the  functional  Ja  {x)  in  (54)  for  the  case  of  this 
CIP.  Let  i7' {Sr)  be  its  Frechet  derivative.  Then  it  was  derived  from  (59)  that 

2 

ll^r  ~  ^r,o(<5)||  ^  ||-^a(5)  11^2(0)  ’ ^ 

In  addition,  the  following  relaxation  property  was  derived  from  Theorem  4.3 

||^r,n+l  ^r,a{5)  ^  ^  ||^r,n  ^r,a(5)  ||  )  h  ^  (b)  1)  j 

as  long  as  7^  ^r,a(S)  and  the  maximal  mesh  step  size  hn+i  on  the  subspace  M„+i  is 
sufficiently  small.  The  estimate  (60)  indicates  that  one  should  rehne  mesh  on  the  rehnement 

step  number  n  +  1  near  those  points  where  the  function  77^(5)  (^r,n)  (x)  .  The  explicit  form 

of  the  function  (^r,n)  (x)  was  obtained  in  a  standard  way  using  the  adjoint  problem 
method.  Thus,  we  arrive  at 

Mesh  Refinement  Recommendation.  Let  (3  G  (0, 1)  he  the  toleranee  number,  which 
is  chosen  numerically.  Refine  the  mesh  in  such  subdomains  of  hi  where 

\K(5)  (^r,n)  (a:)  I  >  (3 1  max  (e^,„)  (a;)|  .  (61) 

In  all  tests  with  adaptivity  below  we  have  used  the  mesh  rehnement  recommendation 
(61). 

5  The  Work  With  Experimental  Data  for  Targets  Placed 
in  Air 

In  this  section  we  describe  our  work  with  backscattering  experimental  data  for  targets  located 
in  the  air.  Results  of  this  section  are  rehected  in  four  publications  [2,  3,  4,  6].  The  main 
challenge  was  a  huge  discrepancy  between  simulated  and  experimental  data,  compare  Figures 
3-a  and  3-b.  Therefore,  one  of  key  difficulties  was  to  invent  a  new  data  preprocessing 
procedure.  The  data  preprocessing  is  inevitably  a  heuristic  procedure.  This  procedure  has 
“moved”  the  pre-processed  data  closer  to  the  simulated  ones  than  the  raw  data.  The  pre- 
processed  data  were  used  as  an  input  for  our  globally  convergent  algorithm.  Therefore,  the 
success  of  our  studies  of  experimental  data  proves  a  signihcant  degree  of  robustness  of  the 
globally  convergent  method  of  this  project. 

We  point  out  that  in  both  preprocessing  of  experimental  data  and  postprocessing  of  calcu¬ 
lated  results  it  is  crucial  to  choose  such  parameters  and  procedures  which  would  work  for  all 
targets  without  any  exceptions.  Or  at  least  one  should  choose  one  set  parameters/procedures 
for  all  dielectric  and  the  second  one  for  all  metallic  targets.  If  such  a  choice  is  possible,  then 
one  can  claim  that  such  a  procedure  is  both  unbiased  and  stable.  This  is  exactly  what  is 
done  below. 
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5.1  Data  acquisition 

Two  main  pieces  of  our  experimental  apparatus  are  Picosecond  Pulse  Generator  and  Tex- 
tronix  Oscilloscope,  see  Figure  1. 


(a)  (b) 

Figure  1:  Two  main  pieces  of  our  experimental  apparatus,  a)  Picosecond  Pulse  Generator 
10070A.  b)  Textronix  Oscilloscope. 


Figure  2:  (a):  A  picture  of  our  experiment  setup;  (b)  Diagram  of  our  setup. 

Our  experimental  configuration  is  shown  in  Figure  2.  Let  E  (x,  t)  =  Ey,  E^)  (x,  t)  be 
the  vector  of  the  electric  held.  Then  it  is  clear  from  Figures  2-a,b  that  u  (x,f)  =  Ey  (x,f) , 
where  u  is  the  function  in  (1).  The  incident  wave  held  is  Ei^c  {x,  t)  =  (0,  5  [x  —  Xq)  6  (f) ,  0) . 
It  was  demonstrated  numerically  in  [29]  that  in  this  case  \Ey  (x,f)|  >>  {E^  {.x,t)  \ ,  \Ez  {x,t)\ 
and  the  propagation  of  the  component  Ey  (x,  t)  is  well  governed  by  equation  (1).  This  and 
especially  the  accuracy  of  our  results  below  justify  the  use  of  equation  (1). 

The  setup  of  our  measurements  included  a  horn  antenna  (transmitter)  hxed  at  a  given 
position  and  a  detector  scanned  in  a  square  of  a  vertical  plane,  which  we  refer  to  as  the 
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measurement  plane.  Consider  the  Cartesian  coordinate  system  Oxyz  as  shown  in  Figure  2b. 
The  scanning  area  was  1  meter  by  1  meter  in  both  horizontal  and  vertical  directions  with 
the  step  size  of  0.02,  starting  at  {x,y)  =  (—0.5,  —0.5),  and  ending  at  {x,y)  =  (0.5,  0.5). 

In  our  mathematical  model  (1),  (2),  we  assume  that  the  source  point  xq  is  in  \ 
n.  However,  due  to  some  technical  difficulties  with  the  mechanical  scanning  system,  the 
horn  antenna  was  not  placed  behind  but  in  front  of  the  measurement  plane  (between  the 
measurement  plane  and  the  targets).  Therefore  a  small  area  in  the  center  of  the  scanning  area 
on  the  measurement  plane  was  shaded  by  the  horn.  The  horn  was  placed  at  the  distance 
of  about  0.2-0.25  from  the  measurement  plane  and  the  distances  from  the  targets  to  the 
measurement  plane  are  about  0.8. 

At  each  position  of  the  detector,  a  number  of  electric  pulses  were  emitted  by  the  horn.  The 
detector  received  two  types  of  signals:  the  direct  signals  from  the  source  and  the  backscatter 
signals  by  the  targets.  The  direct  signals  are  used  for  time  reference  in  data  pre-processing. 
There  were  also  other  unwanted  signals  due  to  scattering  by  some  objects  in  the  room.  To 
reduce  the  instability  of  the  recorded  signals  in  terms  of  magnitude,  the  measurements  were 
repeated  800  times  at  each  detector  position  and  the  recorded  signals  were  averaged.  By 
scanning  the  detector  and  repeating  the  measurements,  we  obtained  essentially  the  same 
signals  as  using  one  incident  signal  and  multiple  detectors  at  the  same  time. 

Pulses  were  generated  by  the  Picosecond  Pulse  Generator  10070A.  The  scattered  signals 
were  measured  by  a  Tektronix  DSA70000  series  real-time  oscilloscope.  The  emitted  pulses 
were  of  300  picoseconds  duration.  The  wavelength  of  the  incident  pulses  was  about  0.03 
m.  The  sampling  rate  (the  step  size  in  time  between  two  consecutive  records  of  captured 
signals)  was  At  =  10  picoseconds.  Each  signal  was  recorded  for  10  nanoseconds. 

5.2  Six  steps  of  data  pre-processing 

One  of  the  biggest  challenges  in  working  with  these  experimental  data  is  the  huge  misfit 
between  these  data  and  the  data  produced  via  computational  simulations.  There  are  several 
causes  of  this  misht  such  as: 

1.  The  instability  of  the  amplitude  of  the  emitted  signals  (incident  waves)  which  causes 
the  instability  of  the  received  signals. 

2.  Unwanted  (parasitic)  scattered  waves  caused  by  the  presence  of  several  existing  objects 
around  our  devices,  see  Figure  2a. 

3.  The  shadow  on  the  measurement  plane  caused  by  the  transmitting  horn  antenna. 

4.  The  difference  between  the  experimental  and  simulated  incident  waves. 

Figure  3  compares  experimental  and  computationally  simulated  signals  at  the  same  de¬ 
tector. 
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A  sample  o1  a  measured  time-dependent  curve 


time  (ns) 


A  sample  of  a  "measured*  time-dependent  curve  for  synthetic  data 


(a) 


(b) 


Figure  3:  Experimental  and  simulated  data  at  the  same  detector,  a)  Experimental  data,  b) 
Computationally  simulated  data. 


Figure  3  compares  time  resolved  experimental  and  simulated  data  at  the  same  detec¬ 
tor.  Computational  simulations  were  performed  for  the  model  target  which  has  the  same 
properties  as  the  one  which  produces  the  signal  on  Figure  3-a). 

Therefore,  the  central  procedure  before  applying  inversion  methods  is  the  data  pre¬ 
processing.  So  that  the  pre-processed  data  would  look  somewhat  similar  with  computa¬ 
tionally  simulated  data.  This  procedure  is  heuristic  and  cannot  be  rigorously  justified.  The 
pre-processed  data  were  used  as  input  for  the  globally  convergent  algorithm.  Our  data  pre¬ 
processing  consists  of  six  steps  described  below.  We  do  not  describe  steps  1-3,  6  in  detail 
here,  since  they  are  straightforward.  However,  steps  4-6  are  completely  new. 

Step  1.  Off-set  correction.  The  acquired  signals  are  usually  shifted  from  the  zero  mean  value. 
This  can  be  corrected  by  subtracting  the  mean  value  from  them. 

Step  2.  Time-zero  correction.  Time-zero  refers  to  the  moment  at  which  the  signal  is  emitted 
from  the  transmitter.  The  recorded  signals  may  be  shifted  in  time.  We  use  the  direct  signals 
from  the  transmitter  to  the  detector  to  correct  the  time-zero. 

Step  3.  Source  shift.  As  mentioned  above,  the  horn  antenna  in  our  experiments  is  placed 
between  the  targets  and  the  measurement  plane.  However,  in  data  calibration,  we  need  to 
simulate  the  data  for  the  case  when  the  measurement  plane  is  between  the  horn  and  the 
targets.  Therefore,  we  artificially  “shift”  the  horn  in  the  positive  z— direction  such  that  it  is 
0.4  m  further  than  the  measurement  plane  from  the  targets.  It  is  done  by  shifting  the  whole 
time-dependent  data  by  a  number  of  samples  which  corresponds  to  the  shifted  distance. 
Step  4-  Extraction  of  scattered  signals.  Apart  from  the  signals  backscatter  by  the  targets, 
our  measured  data  also  contain  various  types  of  signals  as  mentioned  above.  What  we  need, 
however,  is  the  scattered  signals  by  the  targets  only.  The  extraction  of  these  scattered  signals 
for  each  target  is  done  as  follows.  First,  we  exclude  the  direct  signals  and  the  unwanted 
signals,  which  come  earlier  than  the  scattered  signals  by  the  target  (see  Figures  3-a  and  4) 
by  calculating  the  time  of  arrival.  These  unwanted  signals  are  due  to  the  reflection  of  the 
direct  signals  by  the  metallic  structure  behind  the  measurement  plane,  so  we  can  estimate 
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their  times  of  arrival  as  we  know  the  distance  from  the  measurement  plane  to  this  structure. 

We  have  observed  that  the  scattered  signals  by  the  target  are  the  strongest  peaks  of  the 
remaining  ones.  Hence,  after  removing  the  aforementioned  signals,  we  hrst  detect,  for  each 
detector  position,  the  strongest  negative  peak.  Next,  we  consider  the  strongest  negative 
peak.  Suppose  that  the  amplitude  of  the  hrst  negative  peak  prior  the  strongest  is  less  than 
80%  of  the  strongest  one  (see  Figure  3-a)).  Next,  the  scattered  signals  by  the  target  are 
taken  as  7  peaks  (4  negative  peaks  and  3  positive  peaks)  starting  from  that  hrst  negative 
peak  prior  to  the  strongest  one.  However,  if  the  amplitude  of  the  hrst  negative  peak  prior  to 
the  strongest  one  is  greater  or  equals  80%  of  the  strongest  negative  peak,  then  we  start  from 
the  second  negative  peak  prior  to  the  strongest  one.  The  reason  for  choosing  7  peaks  for  the 
scattered  signals  is  due  to  the  fact  that  the  incident  pulses  also  contain  7  strong  peaks.  We 
note  that  having  the  scattered  signals  by  the  target,  we  can  easily  determine  the  distance 
from  its  front  side  to  the  measurement  plane  by  calculating  the  time  of  arrival.  We  use  the 
resulting  signal  for  data  propagation. 

The  next  round  of  extraction  of  scattering  signals  comes  after  data  propagation.  In  this 
case  we  choose  the  largest  negative  peak.  Next,  we  choose  the  closest  negative  peak  to  the 
left  of  it  and  set  to  zero  the  signal  to  the  left  of  this  second  peak.  We  apply  the  Laplace 
transform  (6)  to  the  resulting  curve  to  get  the  boundary  data  for  the  function  w  {x,  s). 

Step  5.  Data  propagation.  After  getting  the  scattered  signals,  the  next  step  of  data  pre¬ 
processing  is  to  propagate  the  data  closer  to  the  targets,  i.e.  to  approximate  the  scattered 
waves  on  a  plane  closer  to  the  targets,  compared  to  the  measurement  plane.  There  are  two 
reasons  for  doing  this.  The  first  one  is  that  since  the  Laplace  transform  decays  exponentially 
in  terms  of  the  time  delay,  which  is  proportional  to  the  distance  from  the  targets  to  the 
measurement  plane,  then  the  amplitude  of  the  data  after  the  Laplace  transform  on  the 
measurement  plane  is  very  small  and  can  be  dominated  by  the  computational  error.  The 
second  reason  is  that  this  propagation  procedure  helps  to  reduce  the  computational  cost 
substantially  as  the  computational  domain  hi  is  reduced.  We  have  also  observed  that  the 
data  propagation  helps  to  reduce  the  noise  in  the  measured  data. 

Step  6.  Data  ealibration.  Finally,  since  the  amplitude  of  the  experimental  incident  and 
scattered  waves  are  usually  quite  different  from  simulations,  we  need  to  bring  the  former  to 
the  same  level  of  the  amplitude  as  the  latter.  This  is  done  using  a  known  target  referred  to 
as  calibrating  object. 

It  is  clear  from  Figure  6  that  the  Laplace  transform  of  the  propagated  data  looks  more 
concentrated  and  less  noisy  than  the  Laplace  transform  of  the  measured  data. 

5.3  Data  propagation 

We  now  describe  in  detail  our  data  propagation  procedure  of  the  above  Step  5  (section  5.2). 
Denote  by  Pm  the  measurement  plane  and  by  Pp  the  propagation  plane,  which  is  closer  to 
the  targets  than  Pm.  Without  a  loss  of  the  generality,  we  denote  by  Pm  =  {z  =  a}  and  Pp  = 
{z  =  0},  where  the  number  a  >  0.  Moreover,  denote  by  F  =  (—0.5,  0.5)  x  (—0.5,  0.5)}  C 
the  scanning  area  of  the  detector  on  the  plane  Pm-  Let  F^  =  {{x,y,a)  G  :  (x,  y)  G  F} 
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Measurement  plane 


Data  propagation  plane 


Figure  4:  Schematic  diagram  of  data  propagation. 

and  Fp  =  0)  G  :  (x,j/)  G  F}.  We  also  denote  by  u^{x,t)  the  scattered  wave. 

Note  that  the  medium  between  and  Pp  is  homogeneous  with  e  =  1  and  the  scattered 
wave  propagates  in  the  direction  from  Pp  to  P^-  The  goal  of  the  data  propagation  is  to 
approximate  ^  from  the  measured  data  q(x,y,t)  := 

irpx(o,oo)  '  ir,nx(o,oo) 

To  do  this,  we  make  use  of  a  time  reversal  method.  Its  idea  is  to  reverse  the  scattered 
wave  in  time  via  solution  of  an  initial  boundary  value  problem  for  the  time-reversed  wave 
function.  We  proceed  as  follows. 

Consider  the  domain  D  :=  {x  =  {x,  y,  z)  E  {x,y)  E  T  ,b  <  z  <  a}  with  6  <  0.  Note 
that  Fp  C  D.  The  reason  for  choosing  this  larger  domain  is  that  we  need  to  assign  boundary 
conditions  at  dD.  Since  short  pulses  are  used  as  incident  waves,  it  is  reasonable  to  assume 
that  the  scattered  wave  in  the  domain  between  Pm  and  Pp  vanishes  along  with  its  time 
derivative  uf  after  some  time  T,  i.e.  (x,  t)  =  0  for  x  E  D,t  >  T.  Therefore,  in  the  following 
we  consider  only  the  finite  time  interval  (0,T).  Denote  t  :=  T  —  t.  Then  the  time-reversed 
wave  function  u^{x,t)  :=  u^{x,t)  satisfies  the  homogeneous  wave  equation.  Moreover,  it 
propagates  in  the  negative  2;  direction,  i.e.  from  Pm  to  Pp.  We  assume  that  the  function 
satisfies  the  absorbing  boundary  condition  at  Ff,  :=  {{x,y,b)  :  {x,y)  E  F}.  On  F;,,  far 
from  our  propagation  plane,  this  boundary  condition  means,  heuristically,  that  we  “send 
back”  the  original  scattered  wave  recorded  at  Pm-  On  the  other  hand,  we  impose  the 
zero  Neumann  boundary  condition  at  the  rest  of  the  boundary  of  D,  except  of  F^.  Denote 
Qt  =  D  X  (0,  T)  and  F3  :=  dD  \  (F^  UF?,).  We  obtain  the  following  problem  for  the  function 


Original  signal  at  detector  ^^1301 


Time  in  nanoseconds 


Extracted  target  signal  at  detector  #1301  after  source  shift 


(a) 


(b) 


Propagated  signal  at  detector  #1254 


Time  in  nanoseconds 


Extracted  target  signal  after  propagation  at  detector  #1 254 


(c)  (d) 

Figure  5:  Extraction  of  scattered  signals,  a)  Original  signal  at  a  detector,  b)  Signal  extracted 
from  a)  after  shifting  the  source,  c)  Propagated  signal.  In  propagation,  extracted  signals  on 
all  detectors  were  used,  c)  Propagated  signal  detected  at  some  point,  d)  Extracted  signal  of 

c) 
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Figure  6:  2D  and  3D  representations  of  the  Laplace  transform  (6)  at  a  certain  value  of 
s  >  0  of  the  data  at  the  measurement  plane  (left  panel)  and  the  propagated  plane  (right 
panel).  It  is  clear  from  this  figure  that  the  Laplace  transform  of  the  propagated  data  looks 
more  concentrated  and  less  noisy  than  the  Laplace  transform  of  the  measured  data. 


t) 


=  Am'’,  (a;,r)  G  Qt, 

{x,  0)  =  u^.  {x,  0)  =  0,  X  E  D, 


u 


'rmxio,T) 


=  9{x,y,T-T), 


{d^vJ'  +  drU 


r\\  ^  Q 

r6x(0,T) 


r3x(o,r) 


(62) 

(63) 

(64) 

(65) 

(66) 


Theorem  1  below  shows  the  stability  of  the  problem  (62)-(66).  This  theorem  is  proven 
via  a  modihcation  of  the  energy  estimates  method  [57],  and  it  can  be  extended  to  more 
general  domains  and  more  general  hyperbolic  operators.  For  brevity  we  are  not  concerned 
here  with  minimal  smoothness  assumptions  and  leave  aside  the  question  of  existence.  We 
conjecture  that  it  can  be  addressed  via  the  technique  of  chapter  4  of  [57]. 


Theorem  1  Assume  that  there  exists  a  solution  E  {Qt)  of  the  problem  (62)-(66). 
Also,  assume  that  the  function  g  E  (T^  x  (0,T))  and  there  exists  such  a  function  F  E 
{Qt)  that 


F  (x,  0)  —  Ft-  (x,  0)  —  0,  {duF  +  Ft-)  |r6x(o,r)—  0,  d^F  |r3x(o,T)—  0, 

F  |r,„x(o,T)=  9  {.X,  t) ,  \\F\\jj2i^q^'j  <  C  ||fi'|li^2('r^x(o,T))  > 

where  C  >  t)  is  a  certain  number.  Then  that  solution  is  unique  and  the  following  stability 
estimate  holds  with  a  constant  Ci  =  Ci  (F,  Qt)  >  0  depending  only  on  the  listed  parameters 

II“’'IIh1(Qt)  -  \\9\\H^(rmX{0,T))  ■ 

By  solving  the  problem  (62)-(66),  we  obtain  an  approximation  of  u'^  {x,t)  and  then 
obtain  an  approximation  of  for  x  EVp.  We  use  the  hnite  difference  method  to  solve  this 
problem. 

5.4  Data  calibration 

Usually  the  experimental  data  have  quite  different  amplitudes  compared  to  the  simulations. 
Figure  7  shows  that  the  minimal  value  of  the  Laplace  transform  of  the  propagated  measured 
data  is  approximately  —2  x  10“^,  whereas  the  minimal  value  for  simulated  data  is  about 
—5  X  10“®.  We  choose  a  number,  which  is  called  calibration  factor,  to  scale  the  measured 
data  to  the  same  scaling  as  in  our  simulations.  To  do  this,  we  make  use  of  the  measured 
data  of  a  single  calibrating  object  whose  location,  shape,  size  and  material  are  known.  The 
word  “single”  is  important  here  to  ensure  that  the  calibration  procedure  is  unbiased,  i.e.  it 
remains  the  same  for  all  targets. 

First,  we  computationally  simulate  the  data  on  Fp  for  the  calibrating  object  by  solving 
the  problem  (68)-(73),  see  section  2.6.  Next,  we  compute  the  Laplace  transform  (6)  of 
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Laplace  transform  of  the  measured  propagated  data,  s  =  10  Laplace  transform  of  the  simulated  data,  s  =  10 


Figure  7:  Laplace  transform  (6)  of  the  scattered  wave  on  the  propagation  plane  Pp.  a)  Mea¬ 
sured  data.  The  maximal  amplitude  of  the  peak  is  2.5  x  10“®.  h)  Simulated  data.  The  maximal 
amplitude  of  the  peak  is  5  x  10“®.  Thus,  the  data  calibration  is  necessary. 

this  computationally  simulated  solution.  Below  we  work  with  s  G  [s,  s] .  Numbers  s,  s  are 
chosen  numerically.  Denote  by  vjh^  (x,  s),  (x,  s)  and  tCgim  (^5  respectively  the  Laplace 

transforms  of  the  total  wave,  the  scattered  wave  and  the  incident  wave  of  the  simulated 
solution  for  the  calibrating  object.  Clearly,  (x,  s)  =  (x,  s)  —  (x,  s).  It  can  be 

proved  that  (x,  s)  <  0.  Figure  7-b)  displays  the  function  (x,  s)  for  x  G  Fp  and 
qualitatively  this  is  a  typical  behavior  for  all  targets.  Let 

dsim,s  =  minw^^^  (x,  s) . 

Fp 

Next,  for  x  G  Fp  let  (x,  s)  be  the  Laplace  transform  of  the  propagated  experimental 
data  for  the  calibrating  object,  see  Figure  7-a).  Denote 

dexp,s  =  min  (x,  s) . 

Fp 

The  number  dsim,s/dexp,s  is  used  as  the  calibration  factor  for  all  targets  at  pseudo-frequency 
s.  That  means,  the  propagated  measured  data  of  all  targets  are  multiplied  by  this  calibration 
factor  before  being  used  in  the  inversion  algorithm. 

We  have  two  types  of  targets:  dielectric  and  metallic  targets.  We  have  observed  that  two 
different  calibration  factors  should  be  used  for  these  two  types  of  targets,  because  the  signals 
from  them  have  different  levels  of  amplitude.  First  of  all,  we  differentiated  these  two  types 
of  targets  by  comparing  the  amplitudes  of  the  recorded  signals.  Indeed,  we  have  consistently 
observed  that  the  maximal  values  of  amplitudes  of  measured  signals  are  at  least  two  times 
larger  for  metallic  targets  than  for  dielectric  ones  on  those  positions  of  detectors  which  are 
most  sensitive  to  the  presence  of  targets.  Next,  we  chose  in  each  type  a  known  object  as 
the  calibrating  object.  In  other  words,  we  should  use  a  dielectric  calibrating  object  for  all 
dielectric  targets  and  another  metal  calibrating  object  for  all  metallic  targets.  As  to  metallic 
targets,  we  have  established  in  [17]  that  one  should  use  the  so-called  appearing  dielectric 
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constant,  whose  values  are 


Er  (metal)  G  [10, 100] .  (67) 

The  value  of  Er  {x)  for  the  dielectric  calibrating  object  was  taken  as  Er  {x)  =  4.28  inside 
that  target  and  e{x)  =  1  outside  of  it.  For  the  metallic  calibrating  object,  as  suggested  by 
(67),  we  took  Er  {x)  =  12  inside  and  e  {x)  =  1  outside  of  it. 

5.5  Dimensionless  variables 

The  spatial  dimensions  in  our  experiment  were  given  in  meters.  Since  the  scanning  step  in 
our  measured  data  was  0.02  m  in  both  x  and  y  directions,  we  chose  the  dimensionless  spatial 
variable  x'  to  be  x'  =  x/l{m).  Next,  to  scale  the  wave  speed  to  be  1  in  the  homogeneous 
medium,  as  in  our  simulations,  we  chose  the  dimensionless  time  t'  by  t'  =  0.3t  where  t 
is  the  time  given  in  nanoseconds  (ns).  The  factor  0.3  is  the  speed  of  light  in  meters  per 
nanosecond  in  the  free  space.  For  the  simplicity  of  notations,  we  use  x  and  t  again  to  denote 
the  dimensionless  variables. 

5.6  Data  simulation 

As  it  is  clear  from  sections  5. 2-5. 4,  we  need  to  computationally  simulate  the  data  via  solving 
the  forward  problem  (1),  (2).  On  the  other  hand,  it  is  clear  that  it  is  impossible  to  solve 
eqnation  (1)  in  the  inhnite  space  R^.  Therefore,  we  have  replaced  R^  with  a  hnite  prism  G 
and  used  absorbing  boundary  conditions  on  those  sides  of  this  prism,  which  are  perpendicular 
to  axis,  and  used  zero  Neumann  bonndary  condition  at  the  rest  of  the  bonndary  of  this 
prism. 

Althongh  a  point  source  is  used  in  the  forward  model  (1),  (2)  for  theoretical  analysis,  we 
make  use  of  an  incident  plane  wave  in  onr  numerical  implementation.  The  prism  G  is 

G  -.=  {x  =  {x,  y,z)  :  Xi  <  X  <  Xu,  Yi  <  y  <  Yu,  Zi  <  z  <  Zu}. 

Denote  by  (9G(,  :=  {z  =  Zi},  dG'^  :=  {z  =  Zu}  and  dGxy  =  dG  \  {dG^^  U  dG^).  An 
incident  plane  wave  of  a  short  time  period  is  excited  at  and  propagates  in  the  negative 
direction.  At  the  plane  we  assnme  that  the  absorbing  boundary  condition  is  satished, 
and  at  dG^y  we  assign  the  homogeneous  Neumann  boundary  condition.  More  precisely,  we 
solve  the  following  problem 


Er{x)Utt{x,t) 

=  Au{x,t),  {x,t)  E  G  X  {0,T), 

(68) 

u{x,  0) 

=  0,  ut{x,0)  =  0,x  E  G, 

(69) 

dyU 

=  f{t),{x,t)  E  dG"}  X  {0,ti), 

(70) 

dyU 

=  -ut,{x,t)  E  dG^  X  {ti,T), 

(71) 

dyU 

=  —ut,{x,t)  E  dG\  X  (0,T), 

(72) 

dyU 

=  t),{x,t)  EdGxyX  {t),T), 

(73) 
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where  v  is  the  outward  unit  normal  vector  of  dG  and  ti  :=  27r/a;  is  the  duration  of  the 
excitation  of  the  incident  plane  wave.  Function  /  is  the  incident  waveform  chosen  by 

fit)  =  sin(a;t),  0  <  t  <  ti  =  ‘In  ju. 

We  have  chosen  ca  =  30.  The  forward  problem  (68)-(73)  is  solved  via  the  hybrid  FEM/FDM 
method  described  in  [1]. 

Because  of  our  specihc  experimental  arrangement,  we  now  specify  the  domain  G.  We 
also  specify  the  computational  domain  hi  where  the  inverse  problem  is  solved.  Given  the 
estimated  distances  from  the  targets  to  the  measurement  plane,  which  were  about  0.8  m, 
we  propagated  the  measured  data  from  the  measurement  plane  Pm  =  {z  =  0.8}  to  the 
propagated  plane  Pp  =  {z  =  0.04}.  This  means  that  the  distance  from  the  front  sides  of 
the  targets  to  the  backscattering  boundary  of  our  inversion  domain  was  about  0.04  m.  The 
reason  for  choosing  this  distance  was  due  to  good  reconstruction  results  we  have  obtained 
for  several  non-blind  targets.  The  domain  G  was  chosen  by 

=  {a;  =  (x,  y,  z)  E  (-0.5,  0.5)  x  (-0.5,  0.5)  x  (-0.1,  0.04)}  .  (74) 

For  solving  the  forward  problem  (68)-(73),  using  the  hybrid  FDM/FEM  method,  we  choose 
the  domain  G  as 

G  =  {x  =  {x,  y,  z)  E  (-0.56,  0.56)  x  (-0.56,  0.56)  x  (-0.16,  0.1)}  .  (75) 

This  domain  G  was  decomposed  into  two  subdomain;  G  =  G  U  (G  \  G).  We  recall  that 
e{x)  =  1  in  G  \  G.  Therefore,  it  is  only  necessary  to  solve  the  inverse  problem  in  G.  In  G 
we  use  a  FEM  mesh  with  tetrahedral  elements,  while  in  G  \  G  we  use  a  FDM  mesh  with  the 
mesh  size  of  0.02  x  0.02  x  0.02  in  Test  1  and  0.01  x  0.01  x  0.01  in  Test  2  below.  The  reason 
for  using  the  FEM  mesh  in  G  is  that  it  is  possible  to  rehne  the  reconstruction  using  adaptive 
mesh  rehnement. 

The  time  interval  (0,T)  in  the  forward  problem  (68)-(73)  was  chosen  equal  to  (0,1.2). 
Since  the  explicit  scheme  in  time  was  used,  the  time  step  size  was  chosen  as  At  =  0.0015 
which  satishes  the  CFL  stability  condition.  The  pseudo  frequencies  Sn  were  chosen  from 
s  =  8  to  s  =  10  with  the  step  size  h  =  0.05.  This  pseudo  frequency  interval  was  chosen 
because  it  gave  good  reconstructions  of  the  non-blind  targets. 

5.7  Complementing  backscattering  data 

Consider  the  boundary  function  iflx,  s)  in  (19).  We  recall  that  only  backscatter  signals  were 
measured  in  our  experiments.  This  means  that  after  data  propagation,  the  function  iflx,  s) 
was  known  only  on  the  side  Fp  =  {x  G  dQ  :  \x\,z  =  0.04}  of  G.  As  in  [16],  we  replace  the 
missing  data  on  the  other  parts  of  dVL  by  the  corresponding  solution  of  the  forward  problem 
in  the  homogeneous  medium.  In  other  words,  in  our  computations,  function  f)  is  given  by 

=  {f^prop{x,s),XE  Tp,SE  (s,s), 

\i^lnif{x,s),x  E  dn\Tp,s  E  {s,s), 
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Figure  8:  Estimation  of  x,y  sizes  of  targets,  a)  Metallic  cylinder,  b)  Wooden  doll  whose 
bottom  is  filled  with  the  sand  and  air  is  in  the  top. 

where  fi’prop  is  the  function  tjj  computed  from  the  propagated  measured  data  at  Fp  and  fiWf 
is  computed  from  the  simulated  solution  of  the  problem  (68)-(73)  with  e{x)  =  1  in  the  wave 
equation  (68). 

Recall  that  our  measured  data  were  collected  with  the  step  size  of  0.02  m  in  a;  and  y 
directions.  To  obtain  the  data  at  the  same  grid  size  as  in  our  computational  domain,  we 
applied  the  linear  interpolation  to  the  Laplace  transform  of  the  propagated  measured  data. 

5.8  Postprocessing 

5.8.1  Estimation  of  the  xy  projection 

We  have  observed  that  the  xy  projection  of  a  target  can  be  roughly  estimated  directly  from 
the  propagated  data.  Let  s  be  the  upper  bound  of  the  pseudo  frequency  interval  on  which 
we  apply  the  globally  convergent  method.  Dehne  F-r  as 

Ft  •  ^propfic  1  y  1  ^propi  s)  ^  0.85  min Up,.op(T,  |/,  s) 

where  Zprop  =  0.04  is  the  value  of  2:  on  the  propagated  plane  Pp  and  Vprop  is  the  function 
u  =  In  (w)  / which  is  constructed  from  the  propagated  measured  data  on  the  propagation 
plane  Fp.  The  truncation  value  0.85  was  chosen  based  on  trial-and-error  tests  on  some  non- 
blind  targets  with  known  sizes.  We  observed  that  Ft  provides  a  good  approximation  for  the 
xy  projection  of  a  target.  Note  that  the  same  truncation  was  applied  to  blind  targets  as 
well.  Hence,  it  is  not  biased. 

5.8.2  Estimation  of  the  z—size  and  the  shape 

The  thing  left  to  estimate  is  the  size  in  the  direction,  i.e.  the  depth  in  the  ;2— direction. 
Along  with  the  above  estimates  of  x,  y— sizes  this  is  supposed  to  provide  an  estimate  of  the 
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shape  of  the  target.  The  z—size  is  more  challenging  than  x,?/— sizes.  At  this  point  of  time 
we  have  a  postprocessing  procedure  which  sometimes  gives  us  a  good  estimate  of  the  2:— size 
and  sometimes  underestimates  it.  We  keep  working  on  rehning  this  procedure.  The  current 
procedure  is  as  follows. 

Let  Srec  {x)  be  the  coefficient  Sr  {x)  reconstructed  by  our  globally  convergent  method.  Let 
Pzo  ■=  {x  =  be  the  plane  where  the  function  Srec  {x)  achieves  its  maximal  value.  Then 
we  compute  the  truncated  coefficient  function  T^ec  {x)  as 


[X 


Srec  (x)  if  Srec  (x,  V,  Zq)  >  7  maX  Erec  {x,  y,  Zq)  , 

1  otherwise, 


(78) 


where  7  G  (0, 1)  is  a  truncation  factor  chosen  such  that  the  area  of  {e{x,y,zo)  >  1}  is  the 
same  as  that  of  Tt,  see  (77)  for  Tt.  Finally,  we  approximate  the  depth  in  the  direction  by 
truncating  Irec  {x)  by  90%  of  its  maximal  value.  This  truncation  value  is  chosen  based  on 
the  trial-and-error  tests  with  non-blind  targets. 


5.9  Targets  in  experiments 

We  have  tested  both  dielectric  and  metallic  targets.  Recall  that  in  metallic  targets  Er  G 
[10,30]  is  the  so-called  appearing  dielectric  constant  of  metals,  see  (67).  In  some  experi¬ 
ments  we  had  only  a  single  target  and  in  others  we  had  two  targets  simultaneously  with  the 
6  centimeters  distance  between  their  surfaces.  We  had  both  non-blind  and  blind  targets. 
“Blind”  means  that  the  computational  group  did  not  have  any  information  about  those  tar¬ 
gets  prior  to  conducting  computations.  Data  for  non-blind  targets  were  used  for  calibration. 
Clearly  blind  cases  are  the  most  challenging  ones.  In  blind  cases  we  have  conducted  compu¬ 
tations  hrst.  And  only  in  a  few  weeks  after  this  we  got  the  true  information  about  targets. 
The  information  was  delivered  by  Mr.  Steven  Kitchin,  an  engineer  who  was  working  part 
time  on  this  grant.  He  has  collected  these  data. 


5.10  Results  of  the  globally  convergent  method 

Tables  1-5  of  this  section  are  quite  informative  ones.  Indeed,  they  provide  the  cumulative 
results  of  our  studies  of  experimental  data. 

5.10.1  Tables 

We  had  total  fourteen  (14)  targets.  Three  (3)  of  them  were  heterogeneous  ones:  (a)  wooden 
doll  with  air  inside,  (b)  the  same  doll  with  a  piece  of  metal  inside  and  (c)  the  same  doll, 
whose  bottom  part  is  hlled  with  sand.  We  believe  that  these  heterogeneous  targets  model 
IDEs,  since  many  of  lEDs  are  mixtures  of  several  substances.  Table  1  provides  a  cumulative 
view  on  these  targets. 
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Table  1:  Targets 


Number  of  targets 

Blind 

Non-blind 

Heterogeneous 

14 

8 

6 

3 

Table  2:  Targets  in  detail.  The  true  information  about  blind  targets  7-14  was  obtained  only 
after  computations  were  conducted.  Blind  targets  12,13,14  are  heterogeneous  ones. 


target  number 

details 

blind 

non-blind 

1 

prism,  oak 

no 

yes 

2 

prism,  pine 

no 

yes 

3 

prism,  pine 

no 

yes 

4 

sphere,  metal 

no 

yes 

5 

cylinder,  metal 

no 

yes 

6 

cylinder,  metal 

no 

yes 

7 

prism,  dielectric 

yes 

no 

8 

prism,  metal 

yes 

no 

9 

two  prisms,  metals 

yes 

no 

10 

prism,  metal 

yes 

no 

11 

two  prisms,  dielectrics 

yes 

no 

12 

wooden  doll,  air  inside 

yes 

no 

13 

wooden  doll,  metal  inside 

yes 

no 

14 

wooden  doll,  sand  inside 

yes 

no 
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Table  2  lists  all  targets  of  out  tests.  Target  3  is  the  same  as  target  2,  but  on  a  different 
distance  from  the  measurement  plane.  Target  6  is  the  same  as  target  5,  but  again  on  a 
different  distance  from  the  measurement  plane. 

Rerfactive  indices  n  =  of  dielectric  targets  were  measured  directly  after  computations 
were  conducted.  We  dehne  the  computed  refractive  index  of  a  target  as 


ncomp  (target)  =  max  J er,comp  (a:) ,  (79) 

Q  V 

where  Sr^comp  [x)  is  the  computed  function  Sr  {x) .  We  have  observed  in  our  computations 
that  the  maximal  value  in  (79)  is  always  reached  within  the  imaged  target,  which  follows 
from  (78).  Table  3  gives  an  information  about  measured  and  computed  refractive  indices  of 
dielectric  targets. 
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Table  3:  Drectly  measured  and  eomputed  refraetive  indices  n  of  targets.  “X/Y”  in  target  11 
means  that  X  corresponds  to  the  upper  target  and  Y  corresponds  to  the  lower  target.  Targets 
7,10,11,12,14  are  blind  ones.  Targets  12,14  are  heterogeneous.  Tests  1,2, 3, 4  differ  from  each 
other  by  different  parameters  of  the  inversion  procedure. 


target  number 

1 

2 

3 

7 

11 

12 

14 

Average  error 

measured 

3.11 

1.84 

1.84 

3.14 

1.84/3.14 

1.89 

3.1 

error 

19% 

18% 

18% 

28% 

18%/28% 

30% 

26% 

23% 

test  1 

1.92 

1.8 

1.81 

1.83 

1.98^-96 

1.86 

1.92 

error 

9% 

3.2% 

1.6% 

16.9% 

7.6%/9.1% 

1.6% 

9.3% 

7.16% 

test  2 

2.01 

2.07 

3.22 

3.21/2.03 

1.83 

3.2 

error 

1.4% 

9.2% 

12.5% 

3.7% 

20.1%/5.4% 

3.2% 

4.8% 

7.5% 

test  3 

2.03 

1.96 

1.65 

3.1 

3.2/3.13 

1.85 

2.05 

error 

3.9% 

6.5% 

11.5% 

1.9% 

19.5%/0.5% 

3.1% 

2.4% 

6.0% 

test  4 

2.02 

2.01 

2.02 

2.03 

2.08/2.06 

1.97 

2.02 

error 

4.4% 

9.2% 

9.8% 

5.4% 

13%/3.9% 

4.2% 

9.6% 

7.43 

THREE  IMPORTANT  OBSERVATIONS  FROM  TABLE  3: 

1.  The  average  error  in  computations  is  at  least  three  (3)  times  less  than  the  average  error 
of  direct  measurements. 

2.  Refractive  indices  of  targets  are  imaged  with  a  good  accuracy,  including  the  most 
challenging  blind  cases  7,11,12,14.  This  is  especially  true  for  cases  of  a  single  target: 
1,2,3,7,12,14. 

3.  The  most  interesting  cases  12,14  of  blind  heterogenous  targets,  which  model  heteroge¬ 
neous  lEDs,  demonstrate  a  good  accuracy:  the  computational  error  is  between  3.1% 
and  9.6%. 
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Table  4:  Computed  appearing  dieleetric  constants  comp  {x)  of  metallic  targets, 

see  (67).  “X/Y”  in  target  9  means  that  X  corresponds  to  the  upper  target  and  Y  corresponds 
to  the  lower  target.  Targets  8,9,13  are  blind  ones.  Target  13  is  heterogeneous.  Tests  1,2, 3, 4 
differ  from  each  other  by  different  parameters  of  the  inversion. 


target  number 

4 

5 

6 

8 

9 

10 

13 

test  1 

14.37 

16.93 

16.45 

25 

12.66/13.1 

13.61 

13.56 

test  2 

15.18 

23.33 

25 

25 

40.53/41.78 

14.13 

14.05 

test  3 

7.59 

10.76 

19.5 

19.5 

11.07/13.1 

8.12 

7.89 

test  4 

15 

15 

15 

15 

13.53/14.06 

15 

14.33 

THREE  IMPORTANT  OBSERVATIONS  FROM  TABLE  4: 

1.  We  can  confidently  image  large  contrasts  in  inclusions.  This  is  hard  to  do  for  locally 
convergent  methods. 

2.  Comparison  of  tables  3  and  4  demonstrates  that  we  can  confidently  distinguish  between 
dielectrics  and  metals. 

3.  The  blind  heterogeneous  target  13  was  a  wooden  doll  with  a  piece  of  metal  inside.  It 
models  the  case  of  an  lED  in  which  an  explosive  is  covered  by  an  innocent  material. 
Nevertheless,  our  calculations  clearly  demonstrate  the  presence  of  a  metal. 

5.11  Some  specifications  for  the  adaptivity 

As  it  was  pointed  out  in  Introduction  and  in  section  4,  we  have  applied  a  two-stage  numerical 
procedure.  The  first  stage  is  the  globally  convergent  numerical  method.  This  stage  has 
provided  accurate  values  of  refractive  indices  n  =  of  dielectric  targets  and  appearing 
dielectric  constants  Sr  of  metallic  targets.  It  is  important  to  point  out  that  the  theory 
guarantees  that  the  first  stage  provides  points  in  a  small  neighborhood  of  the  exact  coefficient. 
Therefore,  the  Tikhonov  functional  does  not  have  local  minima  in  this  neighborhood.  The 
second  stage  is  the  locally  convergent  adaptivity  method.  This  method  uses  the  solution 
obtained  on  the  first  stage  as  the  starting  point  for  iterations.  Adaptivity  refines  images 
via  providing  accurate  shapes  of  targets.  This  is  done  via  sequential  minimization  of  the 
Tikhonov  functional  on  a  sequence  of  locally  refined  meshes  of  finite  elements. 

Data  pre-processing  procedure  for  the  adaptivity  is  a  little  bit  different  from  that  of  the 
first  stage.  This  is  because  the  adaptivity  works  in  time  domain,  whereas  the  globally  con¬ 
vergent  method  works  in  Laplace  transform  domain.  In  addition,  the  image  postprocessing 
procedure  of  the  adaptivity  is  different  from  the  one  of  the  first  stage.  Thus,  we  describe  in 
this  subsection  those  procedures.  We  work  with  propagated  data  and  D,  G,  Tp  are  the  same 
as  in  subsections  6.6,  6.7.  Other  details  can  be  found  in  our  paper  [3]. 
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5.11.1  Data  calibration 

Let  the  function  gexper  {x,t)  ,x  G  Tp,t  G  (0,T)  be  our  propagated  experimental  data.  This 
function  is  given  only  on  grid  points  {xi,  tj) .  We  compute  the  maximal  value  of  this  function, 

9max  max  ffexper  tj)  ■ 


Usually  the  number  g^ax  is  quite  large.  Next,  let  Ugim  {x,  t)  be  the  function  which  is  computed 
via  solving  the  problem  (68)-(73)  where  Sr  (x)  is  the  function  Sr  (x)  computed  by  the  globally 
convergent  method.  Let 

Umax  max  Ugim  t)  ■ 
rpx[o,r] 


Dehne 


Umax 

r  = - 

gmax 


Next,  we  assign 

Uincl  itj)  .  T  ■  5^exper  )  tj ) 

and  use  the  function  gind  {xi,tj)  as  our  data  at  Lp. 


5.11.2  Data  immersing 

The  hnal  step  of  data  preprocessing,  the  so-called  “immersing  procedure”,  which  is  also 
heuristic  and  is  done  as  described  below.  This  procedure  does  two  things: 

•  immerses  the  data  gind  {x,  t)  into  computationally  simulated  ones; 

•  extends  the  data  g^d  {x,t)  from  Tp  to  larger  rectangle  Ti, 

Ti  =  {(x,  y)  G  (-0.56,  0.56)  x  (-0.56,  0.56) ,  x  =  0.04}  . 

To  explain  the  necessity  of  such  data  extension,  we  refer  to  (74),  (75)  and  note  that  the 
domain  12  is  smaller  than  the  domain  G  and  Ti  is  the  orthogonal  projection  of  G  on  the 
plane  {z  =  0.04}.  On  the  other  hand,  we  need  the  data  at  Ti  to  solve  the  so-called  adjoint 
problem  which  appears  in  the  computation  of  the  gradient  of  the  Tikhonov  functional. 

We  dehne  our  immersed  function  Uimmers  {x,t)  for  {x,t)  G  Ti  x  (0,T)  as 

{gind  {x,t) ,  if  a;  G  Tp  and  gind  {x,  t)  >  (3  maxp^  gind  {x,  t) , 

Usim  {x,  t),  if  a;  G  T  and  gind  {x,  t)  <  /3  max^^  gind  {x,  t) ,  (80) 

Usim  (x,  t) ,  if  a:  G  ri^\rp. 

We  choose  the  parameter  (3  =  0.5  in  (80)  in  numerical  experiments  of  section  6.12. 

Figures  9  and  10  show  that,  depending  on  the  parameter  (3  in  (80),  the  data  immersing 
procedure  not  only  allows  to  extend  the  data  from  Tp  to  TiXTp  but  also  make  the  exper¬ 
imental  data  usable  in  our  inverse  algorithm.  Indeed,  we  note  that  the  experimental  data 
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is  measured  at  a  very  high  frequency,  say,  oj  ~  170,  whereas  our  simulations  are  done  at 
O'  =  30  in  order  to  reduce  the  computational  cost.  Therefore,  the  experimental  data  are  not 
compatible  with  the  simulations.  Our  immersing  procedure  helps  to  avoid  solving  the  prob¬ 
lem  at  a  very  high  frequency.  After  this  immersing  procedure  we  solve  the  inverse  problem 
using  the  adaptivity. 


(a)  (3  =  0.2,  f  =  0.3 


60 


Figure  9:  Backscattered  immersed  data  of  the  second  component  E2  of  electric  field  for 
object  number  7  (wooden  doll,  empty  inside)  of  Table  2  for  different  values  of  the  parameter 
(3  in  (80).  Recall  that  the  final  time  is  T  =  1.2. 


5.11.3  Postprocessing  of  results 

Let  Er  (x)  be  the  function  obtained  in  the  adaptive  algorithm.  We  form  the  image  of  the 
dielectric  targets  based  on  the  function  Sr^diei  (x) , 

,  .  f  Srix)  if  Sr  (a;)  >  0.85 maxo^r  (3;) , 

(X)  =  I  j  otherwise. 
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(c)  P  =  0.8,  t  =  0.45 


(d)  (3  =  0.9,  t  =  0.45 


Figure  10:  Backscattered  immersed  data  of  the  second  component  E2  of  electric  field  for 
object  7  (wooden  doll,  empty  inside)  of  Table  2  for  different  values  of  the  parameter  (3  in 
(80).  Recall  that  the  hnal  time  is  T  =  1.2. 


As  to  the  metallic  targets  (i.e.,  the  ones  with  large  computed  maximal  values  of  Sr  (a^)),  we 
use  the  function  Sr^metai  (x) , 


P  _  f  (x)  if  £r  (x)  >  0.3  max^Sr  (x) , 

£r, metal  (xj  |  ^  otherwise. 

We  point  out  that,  in  our  experience,  the  adaptive  algorithm  does  not  change  refractive 
indices  of  dielectric  targets  and  appearing  dielectric  constants  of  metallic  targets. 

5.12  Images  resulting  from  the  two-stage  numerical  procedure 

In  this  subsection  we  present  some  results  of  our  two-stage  numerical  procedure.  The  hrst 
stage  is  the  globally  convergent  method  and  the  second  stage  is  the  adaptivity.  This  pro- 
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Figure  11:  a)  The  image  of  the  target  #1  (wooden  prism,  table  2)  obtained  by  the  globally 
convergent  method,  b)  The  image  of  the  target  #1  obtained  by  the  two-stage  numerical 
procedure,  c)  Zoomed  image  of  b). 


cedure  accurately  reconstructs  all  three  components  of  targets:  refractive  indices,  locations 
and  shapes.  In  the  case  of  metallic  targets  refractive  indices  are  replaced  with  appearing  di¬ 
electric  constants.  While  refractive  indices  and  appearing  dielectric  constants  were  presented 
above  in  Tables  3,4,  here  we  present  samples  of  computed  images  of  shapes  of  targets. 


6  Targets  Buried  in  the  Sand 

We  have  collected  some  experimental  data  for  the  case  of  targets  hidden  in  the  dry  sand. 
Those  are  both  dielectric  and  metallic  targets.  The  data  collection  process  is  continuing. 
Figure  14  displays  our  experimental  arrangement.  Just  as  in  the  experimental  setup  of 
targets  in  the  air,  we  have  The  box  on  this  figure  is  filled  with  the  dry  sand.  We  have 
measured  the  refractive  index  of  sand  n  (sand)  =  2.04.  Hence,  Sr  (sand)  ~  4.16  =  (2.04)^ . 
Note  that  the  front  surface  of  the  sand  box  is  not  flat,  which  complicates  our  inversion 
procedure  even  more.  Targets  are  immersed  in  the  sand  box.  The  distance  between  front 
side  of  the  target  and  the  front  side  of  the  sand  box  varies  between  2  cm  and  10  cm.  Note 
that  we  do  not  know  this  distance  in  advance.  Rather,  we  calculate  it  from  our  data.  To  do 
this,  we  modify  our  data  pre-processing  procedure  of  section  5  for  targets  located  in  the  air. 
The  data  acquisition  process  is  the  same  as  the  one  described  in  subsection  5.1.  The  data 
simulation  is  the  same  as  in  section  5.6. 

First,  we  have  measured  the  signal  from  sand  without  targets,  i.e.  the  reference  signal. 
Next,  we  have  measure  signals  for  the  cases  when  targets  are  immersed  in  the  sand,  see 
Figure  14.  The  major  complication  here  is  that  the  reference  signal  from  the  sand  is  heavily 
mixed  with  the  signal  from  the  case  of  targets  immersed  in  the  sand.  Therefore,  the  data 
preprocessing  procedure  is  again  the  central  one  here  and  it  should  be  different  from  one  of 
the  case  when  targets  are  placed  in  air.  Just  as  in  section  5,  we  had  six  data  preprocessing 
steps: 
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(a)  (b)  (c) 


Figure  12:  a)  The  image  of  the  target  #4  (metallic  sphere,  table  2)  obtained  by  the  globally 
convergent  method,  b)  The  image  of  the  target  #4  obtained  by  the  two-stage  numerical 
procedure,  c)  Zoomed  image  of  b). 


(a)  (b)  (c) 


Figure  13:  a)  The  image  of  the  target  #12  (wooden  doll,  table  2)  obtained  by  the  globally 
convergent  method.  This  is  a  heterogeneous  target  since  air  is  inside  that  doll,  b)  The  image 
of  the  target  #12  obtained  by  the  two-stage  numerical  procedure,  c)  Zoomed  image  of  b). 
One  can  observe  on  c)  that  even  the  “head”  of  the  doll  is  appearing,  although  the  shape  of 
this  doll  is  non  convex. 
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Figure  14:  A  photograph  of  the  experimental  arrangement  of  the  sand  box  measurements. 
Targets  are  immersed  in  this  sand  box.  The  distance  between  the  front  side  of  the  target  and 
the  front  side  of  the  sand  box  varies  between  2  cm  and  1 0  cm.  The  front  surface  of  the  sand 
box  is  not  fiat,  which  complicates  our  inversion  procedure. 

1.  Off-set  correction. 

2.  Time-zero  correction. 

3.  Source  shift. 

4.  Data  propagation. 

5.  Extraction  of  scattered  signals. 

6.  Data  calibration. 

While  steps  7^1,2,3,6  were  basically  the  same  as  ones  in  section  6,  steps  4  and  5  were 
different  and  we  outline  them  here.  We  also  note  that,  unlike  targets  placed  in  air,  where 
different  calibration  factors  were  found  of  metallic  and  dielectric  targets,  we  use  here  only 
one  calibration  factor. 

6.1  A  new  idea  for  data  propagation 

In  section  6  we  have  used  the  time  reversal  method  for  data  propagation,  we  have  observed 
numerically  that  it  works  better  than  a  simpler  Stolt  migration  method  on  computationally 
simulated  data.  However,  on  experimental  data  the  Stolt  migration  method  works  better 
than  the  time  reversal  method.  Basically  the  Stolt  migration  provides  better  quality  data 
than  time  reversal  for  the  case  of  two  targets  present  simultaneously.  The  Stolt  migration  is 
popular  in  Geophysics,  see  [64,  69].  However,  in  the  standard  Stolt  migration  the  wave  at 
the  initial  time  is  calculated  in  the  whole  spatial  domain  of  interest,  whereas  we  calculate 
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the  wave  only  at  a  plane  parallel  to  the  measurement  plane  but  in  the  whole  time  interval. 
This  technique  is  described  as  follows. 

Let  X  =  {x,y,z).  We  assume  that  the  scattered  wave  propagates  in  the  positive  2:- 
direction.  Denote  by  Pm  =  {z  =  b},  b  >  0,  the  measurement  plane  and  by  Pp  =  {z  =  a}, 
with  a  <  b,  the  propagated  plane,  which  is  closer  to  the  target  of  interest  than  Pm-  We 
also  denote  by  u^{x,t)  the  scattered  wave.  Our  objective  here  is  to  determine  the  function 
I/,  t)  :=  u^{x,  y,  a,  t),  given  the  measured  data  /  {x,  y,  t)  :=  u^{x,  y,  b,  t).  We  assume  that 
the  medium  is  homogeneous  in  the  half  space  z  >  a  with  e  =  1.  Therefore,  is  the  solution 
of  the  following  problem: 


=  0,  a;  G  z  >  a,  t  E  {0,  00), 

u'^{x,  0)  =  ul{x,  0)  =  0. 

Consider  the  Fourier  transform 


00  00  00 


u\kx,ky,z,u)  = 


u^{x,y,z,t)e 


—i(ut+xkx-\-yk- 


^^dxdydt. 


0  —00  —00 

It  follows  from  (81)  that  satishes  the  equation: 

Kz  +  -kl-  kl)u‘  =  0,  2  >  a, 

U{k^yky,ZyUj)  ^  ()  (  W)  , 


(81) 

(82) 


(83) 


(84) 

(85) 


where  g{kx,  ky,u)  is  the  Fourier  transform  (83)  of  g{x,y,t).  We  consider  two  cases: 

Case  1:  —  k^  —  ky  <  0.  Keeping  in  mind  that  the  scattered  wave  propagates  in  the 

positive  (^-direction,  the  problem  (84)-(85)  has  the  following  solution 

iilik^,  ky,z,u)  =  g{k^,  ky,u)  exp  (^-(2:  -  a) ^kl  +  k"^-  ,  z  >  a.  (86) 

Case  2:  —  kl  —  ky  >0.  Then  the  solution  -u®  can  be  represented  as 

U2{k^,ky,z,u)  =  g{k^,ky,u)exp  (^-i{z  -kl-  kfj  ,  z  >  a.  (87) 

The  negative  sign  in  the  exponential  term  in  this  formula  is  due  to  the  fact  that  the  scattered 
wave  is  out-going  in  the  positive  ^(-direction. 

Since  the  solution  (86)  is  exponentially  decaying  as  — )■  cxd,  which  represents  the  evanes¬ 
cent  wave,  it  practically  cannot  propagate  to  the  measurement  plane,  which  is  in  the  far  held 
zone.  Hence, 


f{k^,  ky,  u)  =  iilikx,  ky,  b,  u)  =  g{k^,  ky,  u)  exp  (^-i{b  -  a)  -  k^  -  kfj 

Using  the  inverse  Fourier  transform,  we  obtain 
g(x, y,t)=  Iff  f(kz,  ky, 


(88) 


oj^-kl-k^>0 
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Given  the  data  f{x,y,t)  at  the  measurement  plane,  we  compute  /  as  well  as  g{x,y,t)  via 
(88)  using  the  Fast  Fourier  Transform. 

For  each  data  set,  the  propagated  plane  Pp  was  determined  as  follows.  We  hrst  propagated 
the  data  to  the  sand’s  surface.  Using  this  propagated  data,  we  estimated  the  burial  depths 
of  the  targets  (subsection  7.2).  Next,  if  the  burial  depth  of  the  target  closest  to  the  sand 
surface  was  larger  than  4  cm,  we  propagated  the  data  again  from  the  measurement  plane 
up  to  the  plane  Pp,  whose  distance  to  the  front  surface  of  that  target  was  approximately  4 
cm.  Otherwise,  we  used  the  data  propagated  up  to  the  sand’s  surface  for  the  next  step  of 
data  preprocessing.  Note  that  even  we  propagated  the  data  beyond  the  sand’s  surface,  we 
still  saw  the  reflection  from  the  sand’s  surface  in  the  propagated  data  since  we  did  not  take 
into  account  the  presence  of  the  sand  box  in  the  data  propagation,  i.e.,  when  propagating 
the  data  in  the  sand,  we  assumed  that  e  =  1  in  the  sand.  This  reflection  from  the  sand’s 
surface  was  removed  when  the  targets’  signals  were  extracted,  see  subsection  7.2.  Note  that 
the  grid  points  at  Pp  are  the  same  as  the  ones  at  the  measurement  plane  Pm-  Thus,  below 
we  call  “detectors”  the  grid  points  at  the  propagated  plane  Pp. 
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Figure  15:  Result  of  the  data  propagation  for  signals  from  two  targets  buried  inside  the  sand 
box.  The  signals  of  the  two  targets  are  well  separated  from  each  other  as  well  as  from  the 
reflection  from  the  sand’s  surface  after  the  data  propagation. 

A  result  of  the  data  propagation  is  illustrated  in  Figure  15.  The  hgure  shows  a  horizontal 
scan  of  the  sand  box  containing  two  buried  metallic  targets.  The  horizontal  side  denotes  the 
indices  of  the  detector’s  locations  and  the  vertical  side  denotes  time.  Time  increases  from 
the  top  to  the  bottom.  The  propagation  distance,  b  —  a,  is  0.8  m,  which  means  that  the 
propagated  plane  Pp  almost  coincides  with  the  sand’s  surface.  Figure  15(a)  shows  the  original 
data  while  Figure  15(b)  shows  the  data  after  the  propagation.  As  can  be  seen  from  these 
hgures,  the  targets’  signals  in  the  original  data  are  smeared  out.  On  the  other  hand,  they 
are  focused  after  the  data  propagation  making  the  two  targets  more  clearly  distinguished. 
This  is  well  known  for  migration  methods.  Moreover,  we  can  also  see  that  the  reflection  from 
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the  sand’s  surface  is  also  more  visible  after  the  propagation  and  is  well  separated  from  the 
targets’  signals. 

6.1.1  Estimation  of  the  burial  depth 

Since  the  sand’s  surface  reflects  our  microwave  pulses  and  these  reflected  waves  arrive  at  the 
detectors  before  the  ones  reflected  from  the  targets  (see  Figure  15(b)),  the  targets’  signals 
are  dominated  by  that  of  the  sand’s  surface  after  the  Laplace  transform.  In  addition,  the 
measurement  noise  appearing  earlier  than  the  targets’  signals  also  affects  the  latter  after  the 
Laplace  transform  due  to  the  exponential  decay  of  the  kernel,  see  Figure  17(a).  Therefore, 
an  additional  important  data  preprocessing  step  is  needed  in  order  to  pick  up  only  the 
reflection  from  the  targets  and  remove  unwanted  signals  and  noise  coming  earlier  than  the 
targets’  signals.  This  step  is  applied  to  the  propagated  data  as  described  below. 

We  hrst  dehne  some  terms  which  are  used  in  this  subsection.  These  terms  are  related  only 
to  the  propagated  plane  Pp.  The  strongest  detector  in  a  data  set  is  dehned  as  the  detector  at 
which  the  recorded  signal  has  the  largest  amplitude.  A  strong  target  is  either  a  metallic  one 
or  a  nonmetallic  one  whose  dielectric  constant  is  larger  than  that  of  the  sand.  If  the  dielectric 
constant  of  a  target  is  smaller  than  that  of  the  sand,  we  call  it  a  weak  target.  The  strongest 
negative  (positive)  peak  of  a  time-dependent  signal  at  a  certain  detector’s  location  is  the 
negative  (positive)  peak  whose  amplitude  is  larger  than  the  amplitudes  of  other  negative 
(positive)  peaks  of  the  same  signal. 

We  hrst  estimate  the  burial  depth  of  a  target  in  each  data  set.  For  this  purpose,  we  took 
the  strongest  detector.  We  hrst  determined  the  strongest  negative  peak  among  the  hrst  four 
peaks,  starting  from  the  hrst  negative  peak,  see  Figure  16.  This  strongest  negative  peak  is 
considered  as  the  strongest  negative  peak  of  the  sand’s  signal.  After  that,  we  excluded  those 
hrst  four  peaks.  The  reason  for  considering  those  four  peaks  was  due  to  our  observation 
that  those  peaks  should  belong  to  the  rehection  from  the  sand’s  surface.  Moreover,  the  hrst 
two  negative  (so  as  two  positive)  peaks  of  the  incident  wave  were  increasing  in  amplitude. 
After  that,  the  negative  (so  as  positive)  peaks  of  the  incident  wave  decreased  in  amplitude. 
Therefore  any  increase  in  amplitude  of  the  peaks  after  those  hrst  four  peaks  should  be  due 
to  the  rehection  from  the  target.  By  detecting  the  next  negative  (or  positive)  peak  which 
was  stronger  than  the  previous  negative  (positive)  one,  we  located  the  target’s  signal.  Then, 
we  determined  the  strongest  negative  peak  of  the  target’s  signal.  Denoting  by  At  the  time 
delay  between  the  latter  peak  and  the  strongest  negative  peak  of  the  sand’s  signal,  the  burial 
depth  of  the  target  was  approximated  by  n(sand)Af,  where  n(sand)  =  A/£r(sand)  =  2  is  the 
refractive  index  of  the  sand. 

6.2  Extraction  of  target’s  signal:  the  most  difficult  step  of  data 
preprocessing 

After  estimating  the  burial  depth,  we  extracted  the  target’s  signal.  The  extraction  of  signals 
of  targets  in  air  is  quite  simple.  However,  it  is  very  challenging  in  the  case  of  buried  objects. 
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(a)  A  strong  target  (b)  A  weak  target 

Figure  16:  One-dimensional  propagated  signals  at  the  strongest  detectors  of  two  targets: 
one  strong  target  and  one  weak  target.  The  signals  consist  of  the  reflection  from  the  sand’s 
surface  followed  by  the  reflections  from  the  targets. 


especially  for  weak  targets.  Indeed,  when  a  weak  target  is  buried  at  a  shallow  depth,  its 
signal  is  merged  with  the  reflection  from  the  sand’s  surface.  When  it  is  buried  at  a  deep 
depth,  its  signal  is  usually  too  weak  to  be  visible  in  the  data.  Our  experimental  observations 
have  shown  that  if  a  weak  target  is  buried  at  a  depth  of  more  than  5  cm,  then  we  cannot 
detect  it.  In  this  case,  the  target  is  missed. 

As  in  estimating  the  burial  depth,  we  also  worked  with  the  strongest  detector  hrst  and 
excluded  the  hrst  four  peaks.  After  that,  we  selected  the  target’s  signal  as  follows:  (i) 
Suppose  that  either  the  burial  depth  was  larger  than  5  cm,  or  the  strongest  negative  peak  of 
the  target’s  signal  was  larger  than  that  of  the  sand’s  signal  in  amplitude.  Then  we  choose 
as  the  hrst  peak  of  the  target’s  signal  the  strongest  negative  peak  located  after  the  excluded 
ones;  (ii)  otherwise,  the  hrst  peak  of  the  target’s  signal  was  determined  as  the  hrst  positive 
peak  which  was  larger  than  the  previous  positive  one,  provided  that  such  a  peak  exists,  see 
Figure  16.  Since  the  reconstructed  dielectric  constants  of  targets  of  case  (i)  (respectively,  case 
(ii))  was  always  larger  (smaller)  than  that  of  the  sand,  we  also  categorized  a  target  in  case 
(i)  (case  (ii))  as  a  strong  (weak)  target.  For  all  other  detectors,  we  started  from  those  closest 
to  the  strongest  detector  and  on  each  of  them  assigned  as  the  hrst  peak  of  the  target’s  signal 
the  one  which  was  time  wise  closest  to  that  of  the  strongest  detector.  For  strong  targets 
this  one  should  be  a  negative  peak,  while  it  was  a  positive  peak  for  weak  targets.  Next, 
we  continued  similarly  on  all  other  detectors  via  sequentially  choosing  those  peaks  closest 
to  the  one  of  the  previous  detector.  The  reason  for  choosing  a  negative  (positive)  peak  as 
the  hrst  peak  of  the  target’s  signal  for  strong  (weak)  targets  was  due  to  our  observations  in 
numerical  simulations  and  experimental  data  which  have  indicated  that: 

a.  For  a  strong  target,  the  hrst  peak  of  the  target’s  signal  should  be  negative. 

b.  For  a  weak  target,  the  hrst  peak  of  the  target’s  signal  should  be  positive. 
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Moreover,  if  a  strong  target  is  buried  at  a  depth  less  than  5  cm,  then  its  signal  is  stronger 
than  that  of  the  sand  in  amplitude.  If  the  burial  depth  is  more  than  5  cm,  then  its  signal 
might  not  be  stronger  than  that  of  the  sand.  However,  since,  as  we  mentioned  above,  weak 
targets  are  not  visible  at  depths  larger  than  5  cm,  we  consider  all  targets  buried  at  these 
depths  as  strong  targets. 

In  all  above  cases,  the  data  before  the  chosen  hrst  peak  of  the  target’s  signal  were  set  to 
zero.  Hence,  the  Laplace  transform  of  the  preprocessed  data  is  not  affected  by  values  before 
the  hrst  chosen  peak.  We  note  that  such  a  choice  of  starting  peaks  artihcially  immerses 
our  targets  in  air:  because  we  exclude  the  rehection  from  the  sand’s  surface.  Therefore, 
what  we  reconstruct  for  each  target  by  the  globally  convergent  method  is  the  ratio  between 
its  dielectric  constant  (or  the  effective  dielectric  constant  for  metals)  and  that  of  the  sand, 
Er  (target)  /Sr  (sand).  Next,  to  obtain  the  value  of  the  dielectric  constant  of  the  target,  we 
multiply  this  ratio  by  Er  (sand)  =  4. 

Figure  16  shows  one- dimensional  propagated  signals  at  the  strongest  detectors  of  a  strong 
target  and  a  weak  target.  We  indicate  there  the  sand’s  signal  and  the  peaks  of  the  targets. 
These  peaks  were  chosen  as  the  hrst  peaks  of  signals  from  the  targets.  Samples  of  the  Laplace 
transform  of  the  data  before  and  after  the  extraction  of  the  targets’  signals  are  shown  in 
Figure  17,  which  indicate  the  necessity  of  this  preprocessing  step. 

Figure  17  (b)  also  shows  that  the  preprocessed  data  allow  us  to  estimate  locations  of  the 
targets  in  x,y  directions  as  well  as  their  xy-cross  sections,  see  subsection  5.8.1.  These  types 
of  information  help  to  reduce  the  domain  in  which  we  look  for  the  targets.  Indeed,  in  Test 
2  below,  we  took  into  account  these  types  of  information  in  choosing  the  hrst  tail  function. 


(a)  Before  the  extraction  of  the  targets’  signals 
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(b)  After  the  extraction  of  the  targets’  signals 


Figure  17:  The  Laplace  transforms  of  the  data  on  the  propagated  plane  before  (a)  and  after 
(b)  the  extraction  of  the  targets’  signals.  Without  the  extraction,  we  cannot  see  the  targets. 
After  the  extraction,  the  two  targets  show  up  clearly. 
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6.3  Two  tests  and  stopping  criteria 

We  have  analyzed  the  performance  of  the  proposed  algorithm  with  two  different  tests:  Test 
1  and  Test  2.  In  Test  1,  we  made  use  of  the  hrst  tail  function  as  described  in  subsection  3.6, 
with  which  the  global  convergence  is  rigorously  guaranteed.  In  Test  2,  the  estimated  burial 
depth  and  the  xy-cross  section  of  the  target  via  the  data  preprocessing  procedure  were  used 
to  restrict  the  domain  in  which  the  coefficient  8^  was  reconstructed  and  to  choose  the  hrst  tail 
function.  More  precisely,  for  each  target,  let  Xt^min  =  min{a;  G  Ft},  Xt^max  =  max{a;  G  Tr}, 
where  Ft  is  the  estimated  xy  cross  section  of  the  target,  see  subsection  5.8.1.  The  numbers 
yt,min  and  yt,max  are  dehned  similarly.  Then,  we  dehne  the  extended  xy  cross  section  by 

^T,ext  0.03  X  Xf  ^nax  T  0.03,  yt,min  0.03  y  yt,max  T  0.03}. 

Moreover,  denote  by  Ztjront  the  estimated  location  of  the  front  side  of  the  target  in  the  2: 
direction,  given  by  the  burial  depth  estimation.  We  then  dehne  the  following  domain  fir, ext 

^T,ext  •  y  1  G  .  (tC,  I/)  G  0.2  <8  Z  Ziy^QYit  T  0.02}. 

Clearly,  TlT,ext  C  fl.  Moreover,  this  domain  should  contain  the  unknown  target  we  are 
looking  for.  The  last  number  0.02  was  for  compensating  for  possible  error  in  the  estimated 
burial  depth  of  the  target.  Next,  we  chose  the  hrst  tail  function  Vi^i  as  the  function  (41), 
where  the  function  w  {x,s)  was  computed  for  the  coefficient  Er  {x)  :=  Ery  (a^),  where 

Er,0  {x)  =  Er,u,  for  X  G  ^tT,ext,  £r,0  (x)  =  1,  for  X  ^  QT,ext- 

The  upper  bound  Er,u  for  the  function  £^,0  (a^)  was  chosen  as  Er,u  =  25. 

Although  the  convergence  of  the  resulting  algorithm  for  Test  2  has  not  been  rigorously 
proved  yet,  our  numerical  results  show  good  reconstructions,  see  also  [4]  for  accurate  results 
when  targets  are  in  air.  Note  that  we  did  not  use  a  priori  information  about  the  targets. 
Instead,  the  information  used  in  choosing  the  hrst  tail  function  was  derived  from  data  pre¬ 
processing. 

Stopping  criteria: 

Stopping  eriterion  of  Test  1:  The  inner  iterations  with  respect  to  i  are  stopped  at  i  =  ron 
such  that 

Tt n^i  ^  T) ^  ^  imax'j 

where  Dn,i  =  ||fon,i|rp  —  ^r-opl  |L2(rp)-  Here  Tp  is  the  backscattering  side  of  H,  Vprop  is  the  tail 
function  computed  from  the  propagated  data  at  Tp,  and  imax  is  the  maximum  number  of 
inner  iterations.  In  Test  1,  we  have  chosen  imax  =  8. 

The  outer  iteration  with  respect  to  the  pseudo  frequency  s  are  stopped  when  the  error 
function  Dn,i  attains  the  hrst  local  minimum  with  respect  to  n. 

Stopping  eriterion  of  Test  2:  The  inner  iteration  are  stopped  using  the  same  criterion  as 
in  Test  1  but  with  imax  =  5.  The  outer  iteration  is  stopped  when  the  error  function 
i.e.,  the  error  function  at  the  hnal  inner  iteration,  attains  the  hrst  local  minimum. 

We  have  observed  in  our  tests  that  these  stopping  criteria  gave  good  results  for  non  blind 
targets. 
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6.4  Results  of  the  globally  convergent  method  and  some  discus¬ 
sion 

In  Tables  6  and  7  we  summarize  reconstruction  results  of  the  two  tests  for  the  data  sets  listed 
in  Table  5.  Table  6  shows  the  results  for  the  non  metallic  targets.  For  these  targets,  the 
refractive  index  n(target)  =  A/^r (target)  (as  above,  (target)  was  chosen  by  (target)  = 
max^^-^er{x))  is  shown  instead  of  the  dielectric  constant  e  because  n(target)  was  directly 
measured  after  computations  were  performed.  Table  7  shows  the  burial  depths  and  the 
effective  dielectric  constants  of  the  metallic  targets. 

As  described  in  section  6.2,  the  burial  depth  was  estimated  based  on  the  time  delay 
between  the  reflection  by  the  sand’s  surface  and  the  target’s  signal.  Note  that  our  incident 
signal  was  not  really  a  short  impulse.  It  is  therefore  natural  to  expect  some  level  of  error  in 
our  estimates.  Since  we  made  use  of  peaks  of  the  signals  in  estimating  the  depth,  the  error 
we  expect  is  about  the  distance  between  two  consecutive  peaks,  which  is  equal  to  half  of  the 
wavelength  (2  cm). 

From  Tables  6  and  7  we  can  see  that  the  burial  depth  was  accurately  estimated  in  most 
cases,  with  the  errors  not  exceeding  2  cm.  There  are  two  cases  {^4,  7)  in  which  the  errors 
were  about  4  cm.  These  targets  were  buried  at  rather  deep  depths  of  about  the  limiting 
depth  (10  cm)  for  antipersonnel  land  mines.  This  made  the  estimate  less  accurate  because 
of  possible  uncertainty  in  measuring  the  refractive  index  of  the  sand.  Also,  there  might  be 
an  error  in  recording  the  exact  burial  depths  during  the  data  acquisition  for  deeply  buried 
targets. 

The  estimates  of  the  refractive  indices  of  non-metallic  targets  with  refractive  indices  larger 
than  that  of  the  sand  (water  and  wet  wood)  are  quite  accurate  with  the  average  error  of 
about  9.7%  for  Test  1  and  15.2%  for  Test  2.  Note  that  the  error  in  our  direct  measurement 
of  the  refractive  index  of  the  wet  wood  was  10%.  For  water,  we  were  unable  to  directly 
measure  its  refractive  index  at  the  used  quite  high  frequency  of  the  signal,  which  was  about 

7.5  GHz.  Therefore,  we  have  made  a  separate  experiment;  we  have  placed  that  bottle  of 
water  in  air,  measured  the  backscattering  data  and  then  reconstructed  the  refractive  index 
using  the  globally  convergent  method  as  in  [3,  4].  The  result  of  n  =  4.88  matches  well  the 
experimentally  measured  refractive  index  of  4.84  at  high  frequencies  in  Table  3.1  of  [41]. 
Moreover,  by  comparing  our  computed  n  for  water  in  Table  6  with  this  reference  value 
n  =  4.88,  we  can  see  the  consistency  of  our  results. 

Targets  with  smaller  refractive  indices  than  that  of  the  sand  are  of  interest  since  they  are 
models  of  plastic  land  mines  and  lEDs.  We  have  observed  that  we  can  image  these  targets 
only  if  their  burial  depths  do  not  exceed  5  cm.  The  average  error  shown  in  Table  6  for  these 
weak  targets  is  about  21.6%  for  Test  1  and  13%  for  Test  2.  The  average  measurement  error 
of  n  for  weak  targets  was  about  5.4%. 

In  our  experiments,  we  have  missed  some  weak  targets  (not  shown  here),  which  had  more 
than  5  cm  burial  depths.  For  these  weak  targets,  we  have  observed  that  their  signals  were 
blended  by  the  reflection  from  the  sand’s  surface.  Therefore,  we  could  not  detect  any  target’s 
signal  out  of  them.  Note  that,  since  our  current  algorithm  uses  the  Laplace  transform,  it  is 
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Table  5:  Description  of  the  test  data  sets.  Seven  of  them  consist  of  two  targets  each  (5,  6, 
16,  17,  20,  23,  25).  Two  targets  can  he  considered  as  heterogeneous  (11,  12). 


Object 

# 

Blind/ 

Non-blind 

Description  of  target 

Material 

1 

Non-blind 

A  metallic  cylinder 

Metal 

2 

Non-blind 

A  metallic  ball 

Metal 

3 

Non-blind 

A  bottle  filled  with  clear  water 

Water 

4 

Non-blind 

A  wet  wooden  block 

Wet  wood 

5 

Non-blind 

Two  metallic  blocks  at  6  cm  separation 

Metal/Metal 

6 

Non-blind 

A  metallic  cylinder  and  a  teflon  bar 

Metal/Teflon 

7 

Non-blind 

A  metallic  block 

Metal 

8 

Non-blind 

An  empty  bottle 

Air 

9 

Non-blind 

A  bottle  filled  with  teflon  bars 

Teflon 

10 

Blind 

A  ceramic  mug 

Ceramic 

11 

Blind 

A  wooden  doll  filled  with  metallic  screws 
(heterogeneous,  diffuse  scattering) 

Wood/Metal 

12 

Blind 

A  geode  (heterogeneous): 

two  spherical  layers  and  air  inside 

Rock 

13 

Blind 

A  piece  of  rock 

Rock 

14 

Blind 

A  plastic  bottle  filled  with  coffee  grounds 

Coffee  grounds 

15 

Blind 

A  ceramic  mug 

Ceramic 

16 

Blind 

A  cylinder  and  a  block  at  3  cm  separation 

Metal/Metal 

17 

Blind 

An  aluminum  can  and  a  block 

Metal/Metal 

18 

Blind 

A  wooden  doll  with  a  metallic  block  inside 
(heterogeneous) 

Wood/Metal 

19 

Blind 

A  bottle  of  water 

Water 

20 

Blind 

A  metallic  block  and  a  rock 

Metal/rock 

21 

Blind 

A  steel  mug 

Metal 

22 

Blind 

A  wet  wooden  block 

Wet  wood 

23 

Blind 

A  wet  wooden  block  and  an  empty  bottle 

Wet  wood/ air 

24 

Blind 

A  wet  wooden  block 

Wet  wood 

25 

Non-blind 

Two  metallic  blocks  at  1  cm  separation 

Metal/Metal 
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Table  6:  Result  of  the  globally  convergent  algorithm:  the  refractive  indices  n  =  i/e  and 
the  burial  depths  of  non-metallic  targets.  Object  #11  is  a  heterogeneous  target  with  diffuse 
scattering,  see  below.  Object  #12  is  a  heterogeneous  one  with  outer  and  inner  layers,  the 
computed  n  is  compared  with  the  average  measured  n  =  1.28.  Object  #23  consists  of  two 
targets:  wet  wood  and  empty  bottle  filled  with  air.  “Comp.”  stands  for  “Computed”.  The 
average  error  of  strong  targets  is  8.5%  for  Test  1  and  14.7%  for  Test  2.  The  average  error 
of  weak  targets  is  21.6%  for  Test  1  and  13%  for  Test  2. 


Object 

# 

Material 

Comp. 

depth 

Exact 

depth 

Comp, 
n,  Test  1 

Comp, 
n,  Test  2 

Measured 

n 

3 

Water 

3.6  cm 

4.0  cm 

4.7 

4.9 

4.88 

4 

Wet  wood 

5.5  cm 

9.8  cm 

4.4 

4.5 

4.02 

8 

Air 

2.8  cm 

3.0  cm 

1.0 

0.98 

1.0 

9 

Teflon 

2.9  cm 

2.5  cm 

1.0 

1.18 

1.0 

10 

Ceramic 

4.0  cm 

5.0  cm 

1.0 

1.23 

1.39 

11 

Wood  with 
metal  screws 

4.6  cm 

4.0  cm 

1.0 

1.46 

1.89  (wood) 

N/A:  diffuse  scattering 

12 

Geode 
(two  layers) 

3.1  cm 

2.5  cm 

1.0 

1.52 

1.31  (outer) 

1.25  (inner) 

1.28  (average) 

13 

Rock 

2.0  cm 

2.3  cm 

1.0 

1.34 

1.34 

14 

Coffee  grounds 

2.0  cm 

2.5  cm 

1.0 

1.46 

1.11 

15 

Ceramic 

2.6  cm 

2.5  cm 

1.0 

1.51 

1.39 

19 

Water 

7.5  cm 

9.5  cm 

4.5 

5.2 

4.88 

22 

Wet  wood 

2.9  cm 

3.0  cm 

4.8 

5.3 

4.02 

23 

Wet  wood 

5.7  cm 

7.5  cm 

4.0 

4.1 

4.02 

Empty  bottle  (air) 

missed 

7.5  cm 

missed 

missed 

1.0 

24 

Wet  wood 

5.1  cm 

6.8  cm 

3.67 

3.0 

4.02 
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Table  7:  Result  of  the  globally  convergent  algorithm:  the  estimated  effective  dielectric  con¬ 
stants  and  the  burial  depths  of  metallic  targets.  Object  #18  is  a  heterogeneous  one:  a  wooden 
dool  with  a  metallic  block  inside.  Object  #20  consists  of  two  targets:  a  metallic  block  and 
a  rock.  Measured  n{rock)  =  1.34.  Objects  #5,  16,  17  consist  of  two  metallic  targets.  Ob¬ 
ject  #25  consists  of  two  metallic  targets  with  1  cm  distance  between  their  surfaces:  super 
resolution,  see  Figure  B.f. 


Object 

# 

Material 

Computed 

depth 

Exact 

depth 

Computed  e 
Test  1 

Computed  e 
Test  2 

1 

Metal 

2.9  cm 

4.0  cm 

29.9 

46.4 

2 

Metal 

2.9  cm 

3.0  cm 

24.5 

31.0 

5 

Metal 

3.0  cm 

3.0  cm 

23.4 

32.4 

Metal 

3.6  cm 

3.0  cm 

30.5 

41.2 

6 

Metal 

Teflon 

2.8  cm 

missed 

8.5  cm 
8.5  cm 

27.8 

37.5 

7 

Metal 

9.9  cm 

14.0  cm 

47.4 

65.8 

16 

Metal 

2.5  cm 

4.5  cm 

19.9 

24 

Metal 

3.7  cm 

4.5  cm 

33.7 

47.5 

17 

Metal 

2.0  cm 

3.8  cm 

30.0 

51.1 

Metal 

2.7  cm 

3.8  cm 

54.8 

93.5 

18 

Wood,  metal  block  inside 

7.1  cm 

8.5  cm 

18.3 

19.9 

20 

Metal 

Rock 

6.8  cm 
missed 

8.5  cm 
8.5  cm 

30.0 

48.1 

21 

Metal 

5.1  cm 

7.5  cm 

23.1 

28.2 

25 

Metal 

3.8  cm 

4.0  cm 

70.0 

99.8 

Metal 

4.0  cm 

4.0  cm 

40.8 

56.5 
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applicable  only  when  we  can  detect  and  extract  targets’  signals  and  remove  the  reflection 
by  the  sand’s  surface  as  well  as  noise  at  earlier  times.  Otherwise,  they  will  dominate  the 
targets’  signals  after  the  Laplace  transform,  see  Figure  17.  Thus,  these  missed  cases  were 
not  due  to  the  inversion  algorithm. 

The  signals  of  the  metallic  targets  were  strong  compared  to  the  sand’s  signal.  Therefore 
they  were  quite  easy  to  detect.  We  recall  that  metallic  targets  can  be  approximated  by 
dielectric  ones  with  large  appearing  dielectric  constants.  The  estimated  effective  dielectric 
constants  of  our  metallic  test  targets  are  between  20  and  100. 

In  our  tests,  there  were  four  cases  (^^5,  16,  17,  25  in  Table  7)  in  which  there  were  two 
metallic  targets  simultaneously.  In  each  case  we  have  accurately  imaged  both  targets. 

There  were  also  three  data  sets  (7^23  in  Table  6  and  7^6,  20  in  Table  7)  in  which  one 
strong  target  and  one  weak  target  were  buried  simultaneously.  In  all  these  three  cases  we 
have  accurately  imaged  the  stronger  target.  However,  we  missed  the  weaker  ones.  The 
reason  why  we  missed  the  weak  targets  was  due  to  the  fact  that  their  burial  depths  were 
larger  than  5  cm,  which  is  our  limiting  depth  for  weak  targets. 

Also  of  interest  are  three  cases  of  heterogeneous  targets  (7^1 1,  12  in  Table  6  and  ^^18  in 
Table  7),  since  explosive  devices  are  heterogeneous  sometimes.  We  successfully  estimated 
the  average  refractive  index  of  the  geode,  which  consists  of  two  different  layers,  in  data  set 
7^12.  For  the  wooden  doll  containing  a  metallic  block  inside  in  data  set  7^18  the  computed 
dielectric  constant  is  larger  than  that  of  the  wood  but  smaller  than  other  metallic  targets.  It 
is  smaller  because  the  wood  covers  the  metal.  Target  7^11  was  a  wooden  doll  with  randomly 
distributed  metal  screws  inside.  In  this  case  we  observed  a  weak  signal,  rather  than  a  strong 
one  from  the  metal.  In  fact,  we  observed  a  well  known  phenomenon  of  dijfuse  seattering, 
which  was  described  in  [62].  This  can  be  explained  since  the  metal  screws  were  randomly 
oriented  and  represent  a  conducting  very  rough  surface  to  the  incident  microwave  pulse. 
Multiple  scattering  combined  with  the  penetration  of  the  microwaves  into  the  gaps  between 
the  screws  strongly  attenuates  the  incident  wave  and  little  scatters  back  to  contribute  to  a 
measurable  signal. 

6.5  Some  images 

In  this  subsection  we  present  samples  of  images  of  buried  objects  obtained  via  the  application 
of  our  two-stage  numerical  procedure.  The  first  stage  is  the  globally  convergent  method  and 
the  second  stage  is  the  adaptivity.  Just  as  above  (section  5.12),  this  procedure  accurately 
reconstructs  all  three  components  of  targets:  refractive  indices,  locations  and  shapes.  The 
choice  of  parameters  and  calibration  procedure  here  are  the  same  as  in  subsection  5.11.  Left 
panels  of  Figures  18-21  are  images  obtained  on  the  first  stage  and  right  panels  are  images 
obtained  on  the  second  stage.  For  a  better  visualization,  we  zoom  images  via  showing 
reconstruction  results  in  the  box  0.4x0.4x0.24. 
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Figure  18:  Reconstruction  of  a  metallic  ball  (object  in  Table  5).  a)  Result  of  the  globally 
convergent  method  (stage  1).  b)  Result  of  the  adaptivity  (stage  2).  One  can  observe  that  the 
shape  is  reconstructed  rather  accurately  on  b). 


6.6  Super  resolution  the  first  experimental  and  numerical  obser¬ 
vations 

In  the  data  set  7^25,  the  two  metallic  blocks  were  at  1  cm  distance,  see  Figure  21.  On 
the  other  hand,  the  wavelength  of  our  device  is  about  4  cm.  Thus,  the  super  resolution  is 
achieved,  which  is  of  about  A/4,  where  A  is  the  wavelength  of  our  incident  wave.  This  is  an 
unexpected  surprise.  From  a  purely  angular  spectrum  argument,  the  spread  of  backscatter 
angles  for  a  hxed  frequency  would  suggest  a  resolution  of  half  a  wavelength.  There  has  been 
previous  evidence  reported  that  using  a  nonlinear  inverse  scattering  algorithm  for  which 
strong  or  multiple  scattering  occurs,  that  some  degree  of  super  resolution  (i.e.,  beyond  the 
ideal  diffraction  limit  of  half  a  wavelength)  can  occur,  see,  e.g.,  [63].  Still,  this  was  done  in 
[63]  only  for  computationally  simulated  data. 

However,  we  are  the  first  ones  who  has  observed  this  phenomenon  experimentally  and, 
at  the  same  time,  reconstructed  the  image  numerically. 
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(a) 


(b) 


Figure  19:  Reconstruction  of  a  plastic  bottle  filled  with  a  clean  water  (object  in  Table 
5).  a)  Result  of  the  globally  convergent  method  (stage  1).  b)  Result  of  the  adaptivity  (stage 
2).  Since  the  target  was  guite  high  (about  18  cm  tall),  then  the  incident  signal  was  rather 
weak  at  the  top  and  bottom  of  that  bottle.  This  made  it  hard  to  reconstruct  the  shape  more 
accurately.  Still,  one  can  observe  a  stretch  in  the  vertical  direction. 
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(a)  (b) 

Figure  20:  Reconstruction  of  a  ceramic  mug  (object  #10  in  Table  5).  a)  a)  Result  of  the 
globally  convergent  method  (stage  1).  b)  Result  of  the  adaptivity  (stage  2).  One  can  observe 
that  the  shape  is  reconstructed  rather  accurately  on  b). 


7  Joint  Research  With  Engineers  of  US  Army  Research 
Laboratory  and  Technology  Transfer 

Publications  [17,  18,  19]  are  joint  ones  with  Doctors  Anders  Sullivan  and  Lam  Nguyen, 
engineers  of  RF  Signal  Processing  and  Modeling  Branch  of  Sensors  and  Electron  Devices 
Directorate  of  US  Army  Research  Laboratory  (ARL).  We  have  obtained  the  experimental 
data  collected  in  US  Army  Research  Laboratory  (ARL).  The  data  were  kindly  delivered  to 
us  by  Drs.  Anders  Sullivan  and  Lam  Nguyen.  These  data  were  collected  by  the  forward 
looking  radar  of  ARL  [59].  The  goal  of  this  radar  is  to  detect  and  possibly  identify  shallow 
explosive-like  targets,  which  look  like  explosives,  including  lEDs.  Shallow  targets  can  either 
lie  on  the  ground  or  a  few  centimeters  below  the  ground.  Our  goal  was  to  image  ratios  R, 

^  Er  (target) 

Er  (background)  ’ 

using  our  globally  convergent  numerical  method.  If  the  dielectric  constant  of  the  background 
is  known,  then  the  dielectric  constant  of  the  target  can  be  calculated  from  (89),  as  long  as 
R  is  known.  These  were  completely  blind  data.  In  other  words,  we  had  no  idea  what  kind 
of  targets  were  in.  The  only  thing  we  knew  was  that  a  target  of  interest  was  located  either 
above  the  ground  or  buried  in  the  ground.  In  addition,  there  were  many  other  uncertainties 
in  these  data.  A  signihcant  difficulty  was  that  targets  were  hidden  in  a  cluttered  background. 
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(a)  (b) 

Figure  21:  Super  resolution.  Reconstruction  of  two  metallic  targets  (object  ff25  in  Table 
5).  a)  Result  of  the  globally  convergent  method  (stage  1).  b)  Result  of  the  adaptivity  (stage  2). 
One  can  observe  that  the  shapes  are  reconstructed  rather  accurately  on  b).  The  wavelength 
of  our  signal  was  A  =4  cm.  On  the  other  hand,  the  distance  between  surfaces  of  these  two 
targets  was  1  cm.  Therefore,  this  is  a  very  rare  experimental  and  numerical  observation  of 
the  super  resolution  phenomenon. 
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In  the  paper  [19]  we  have  compared  performances  of  two  numerical  methods  on  the  set 
of  experimental  data  or  ARL; 

1.  The  ID  version  of  the  globally  convergent  method  of  this  project. 

2.  The  classical  ID  method  of  the  Krein  equation,  which  is  known  since  1954  [51]. 

Detailed  explanations  of  the  Krein  equation  method  can  be  found  in  [36,  46].  It  was 
demonstrated  in  [19]  that  while  our  method  computes  the  value  of  R  in  (89)  pretty  accu¬ 
rately,  the  Krein  equation  method  completely  fails  on  the  ARL  data.  Furthermore,  the  poor 
performance  of  the  Krein  equation  method  was  explained  analytically,  see  subsection  7.4. 

7.1  Technology  transfer  to  ARL 

We  have  created  an  easy-to-use  software  which  works  with  the  above  ARL  data.  Next,  we 
have  transferred  this  software  to  Drs.  Sullivan  and  Nguyen. 

7.2  Results  of  [17,  18] 

To  image  ratios  (89),  we  have  applied  the  ID  version  of  our  globally  convergent  numerical 
method  of  this  project.  Next,  Drs.  Sullivan  and  Nguyen  have  compared  our  results  with 
those  in  tables  [65,  66]  (dielectric  constants  were  not  measured  in  experiments).  After  these 
comparisons,  both  mathematical  and  engineering  teams  were  quite  happy  that  computational 
results  were  within  tabulated  limits,  see  Table  8.  This  points  towards  both  the  adequacy 
and  the  robustness  of  our  algorithm. 

The  recovered  dielectric  constant  by  itself  is  not  a  sufficient  information  to  distinguish  one 
target  from  another.  The  purpose  of  estimating  the  dielectric  constant  is  to  provide  one  extra 
piece  of  information  about  the  target.  Indeed,  up  to  this  point,  most  of  the  radar  community 
relies  solely  on  the  intensity  of  the  radar  image  for  doing  detection  and  discrimination  of 
targets.  It  is  hoped  therefore  that  when  the  intensity  information  is  coupled  with  the  new 
dielectric  information,  algorithms  could  then  be  designed  that  will  ultimately  provide  better 
performance  in  terms  of  probability  of  detection  and  false  alarm  rate.  As  it  is  clear  from 
Table  8,  some  targets  will  have  dielectric  values  that  tend  to  group  together,  but  even  that 
is  a  useful  information.  For  example,  if  the  estimated  dielectric  value  is  consistent  with  a 
plastic  land  mine,  then  this  would  be  another  clue  to  uncovering  the  target. 

Let  a;  G  M  be  the  free  variable  in  the  ID  case  of  this  study.  In  our  mathematical  model 
the  ratio  R  in  (89)  depends  on  x,  R  =  R  (x) ,  since  we  treat  R  (x)  as  the  unknown  coefficient 
Er  (x)  in  the  ID  analog  of  equation  (1).  Let  Rcomp  (x)  be  function  R(x)  computed  by  our 
globally  convergent  method.  Then  we  define  the  target/background  contrast  R  as 

_  f  Rcomp  (x)  if  Rcomp  (x)  ^  1, 

\  min  Rcomp  (x)  if  Rcomp  (x)  <  1. 
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Table  8:  Blindly  computed  by  the  globally  convergent  method  dielectric  constants  for  five  avail¬ 
able  data  sets  of  targets  immersed  in  cluttered  backgrounds  [17,  18],  “tabulated  Sr  {target)” 
means  either  Er  from  (89)  for  metallic  targets,  or  Sr  from  tables  [65,  66]  for  dielectrics,  or 
Sr  from  [39]  for  bush.  “1”  in  the  second  column  means  that  the  target  was  in  the  air,  and 
“[3,5]”  means  that  it  was  buried  in  the  dry  sand  whose  Er  {dry  sand)  G  [3,5]  [65]. 


target 

R 

Er  (background) 

computed  Er  (target) 

tabulated  Er  (target) 

metal  cylinder 

4.3 

|3,5] 

[12.9,21.4] 

[10,30] 

metal  box 

3.8 

|3,5] 

[11.4,19] 

[10,30] 

wood  stake 

3.8 

1 

3.8 

[2,6] 

bush  (clutter) 

6.4 

1 

6.4 

[3,20] 

plastic  cylinder 

0.28 

|3,5] 

[0.84,1.4] 

1.2 

CONCLUSION  FROM  TABLE  8; 

Dielectric  constants  were  not  measured  in  these  experiments.  Therefore,  the  crucial 
conclusion  from  this  table  is  that  all  computed  values  of  dielectric  constants  of  targets 

are  within  tabulated  limits.  This  was  done  for  the  most  challenging  case  of  blind 
experimental  data  for  targets  immersed  in  cluttered  backgrounds. 

7.3  Data  pre-processing 

To  obtain  results  listed  in  Table  8,  a  new  data  preprocessing  procedure  was  applied  in  [18,  17]. 
This  was  necessary,  because  the  experimental  data  look  very  much  different  from  simulated 
ones,  just  as  in  sections  5,6.  So  that  the  pre-processed  data  would  look  somewhat  similar 
with  computationally  simulated  data.  The  pre-processed  data  were  used  as  input  for  the 
globally  convergent  algorithm. 

The  data  processing  consisted  of  the  following  two  stages: 

Stage  1.  Selection  of  a  single  peak  out  of  the  entire  experimental  curve  by  the  following 
rule:  this  should  be  the  earliest  peak  of  the  largest  amplitude 

out  of  /  peaks  for  a  target  buried  in  the  ground, 

[  all  downwards  looking  peaks  for  a  target  above  the  ground. 

Stage  2.  The  choice  of  the  calibration  factor  CF.  This  was  necessary  because  the 
amplitudes  of  the  experimental  data  were  of  the  order  of  10^  whereas  amplitudes  of  com¬ 
putationally  simulated  data  were  less  than  1.  We  have  chosen  CF  in  such  a  way  that  the 
Laplace  transform  of  the  simulated  data  was  rather  close  to  the  Laplace  transform  of  the 
those  selected  peaks  of  experimental  data,  when  those  peaks  were  multiplied  by  CF.  Thus, 
we  got  CF  =  10“^  and  our  pre-processed  data  were 

pre-processed  data=10“^  ■  (data  pre-processed  on  the  first  stage) .  (90) 

Remark.  The  contrast  i?  of  a  target  heavily  depends  on  the  calibration  factor  CF. 
Therefore,  it  is  very  important  to  have  a  single  calibration  factor  for  all  targets.  Otherwise, 
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Stage  2  would  be  biased.  We  point  out  that  we  did  not  know  R  in  advance  because  of  the 
blind  study  case. 

Figure  22  shows  original  and  pre-processed  data. 


(a)  (b) 

Figure  22:  Original  and  pre-processed  Army  Research  Laboratory  data,  a)  A  sample  of  the 
original  time  resolved  data  data,  b)  Pre-processed  data  for  all  five  available  targets  (super¬ 
imposed).  The  data  pre-processing  was  done  as  described  in  Stages  1  and  2. 


7.4  Comparison  with  the  classical  Krein  equation  method  [19] 

The  classical  Krein  integral  equation  is 

Z 

2/(+0)p(2:,t) -b  j  f'{t  -  \)p{z,\)d\  =  -1.  (91) 

—  Z 

Here  z  is  a  free  variable,  which  is  connected  with  the  original  variable  x  via  {x) ,  where 
2:  {x)  is  a  one-to-one  function.  Next,  /  (t)  is  the  data  function  which  is  extended  as  an  odd 
function  from  {t  >  0}  into  {t  <  0} .  The  function  p{z,  t)  depends  on  2:  as  on  a  parameter. 
So  that  equation  (91)  is  solved  for  each  z  >  0.  As  soon  as  the  function  p{z,t)  is  fonnd,  the 
nnknown  coefficient  R  {x)  can  be  easily  reconstrncted. 

We  took  the  preprocessed  data  of  the  above  Stage  1  and  have  tried  to  find  a  calibration 
factor  CFK  for  them,  in  snch  a  way  that,  similarly  with  (90),  we  wonld  have  for  all  five 
targets 

pre-processed  data  =  f  (t)  =  CFK  ■  (data  pre-processed  on  the  hrst  stage) ,  (92) 

see  Remark  in  snbsection  7.3.  To  find  the  nnmber  CFK,  we  have  selected  any  of  above  hve 
targets.  We  call  it  calibration  target.  Then  have  fonnd  a  “temporary”  CFiF(target)  in  snch 
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Table  9:  Computed  values  of  Rk  for  Case  1.  Wood  stake  is  the  ealibration  target.  Bold 
faced  are  values  of  Rk  whieh  lead  to  values  of  Sr  (target)  outside  of  tabulated  intervals. 


target 

Rk 

Sr  (bckgr) 

computed  Sr  (target) 

tabulated  Sr  (target) 

metal  cylinder 

3.25 

|3,5] 

[9.75,16.25] 

[10,30] 

metal  box 

3.97 

|3,5] 

[11.91,19.85] 

[10,30] 

bush 

42.5 

1 

42.5 

[3,20] 

plastic  cylinder 

0.6 

|3,5] 

[1.8,3] 

1.2 

Table  10:  Computed  values  of  Rk  for  Case  2.  Plastie  cylinder  is  the  ealibration  target.  Bold 
faced  are  values  of  Rk  whieh  lead  to  values  of  Sr  (target)  outside  of  tabulated  intervals. 


target 

Rk 

Er  (bckgr) 

computed  Er  (target) 

tabulated  Er  (target) 

metal  cylinder 

9.8 

|3,5| 

[29.4,49] 

[10,30] 

metal  box 

13.9 

|3,5| 

[41.7,69.5] 

[10,30] 

wood  stake 

17 

1 

17 

12.6] 

bush 

»100 

1 

»100 

[3,20] 

a  way  that  the  value  of  i?  ;=  Rk  for  this  target  computed  via  equation  (91)  would  be  the 
same  as  the  one  listed  in  Table  9.  Next,  we  have  multiplied  the  data  for  four  other  targets 
by  the  same  number  CFK(taTget)  and  have  solved  Krein  equation  (91)  for  so  pre-processed 
data.  Next,  we  have  done  the  same  for  another  target  and  repeated  this  five  times. 

We  have  found  that  it  is  impossible  to  choose  a  proper  number  CFif  (target)  in  such  a 
way  that  resulting  values  of  Rk  would  be  within  tabulated  limits  for  all  five  targets.  More 
precisely,  for  each  out  of  five  possible  CFi^(target)  there  existed  at  least  one  target  for 
which  the  value  of  Rk  was  outside  the  tabulated  limit.  Below  are  three  cases  illustrating 
this  conclusion.  In  the  case  when  metal  box  was  chosen  as  the  calibration  target,  our 
calculations  led  to  the  same  results  as  those  in  Table  6.  This  is  because  R  (wood  stake)  = 
i?  (metal  box)  =  3.8.  Also,  i?  (metal  cylinder)  =  4.3  ~  3.8  =  i?  (metal  box) .  Hence,  the 
calculation  for  the  case  when  the  metal  cylinder  was  chosen  as  the  calibration  target  again 
led  to  results,  which  are  similar  with  those  of  Table  9. 

Case  1.  Wood  stake  is  the  calibration  target.  We  have  obtained  CFK  (wood  stake)  = 

10-5. 

Case  2.  Plastic  cylinder  is  the  calibration  target.  We  have  obtained  C'FiP(plastic  cylinder) 
1.8  ■  10-^ 

Case  3.  Bush  is  the  calibration  target.  We  have  obtained  C'FiP(bush)  =  0.6  ■  10“®. 

To  see  how  the  changes  in  the  calibration  factor  affect  the  reconstructed  values  of  R  and 
Rk,  we  present  Figure  23.  The  visual  analysis  of  these  curves  shows  that  the  increase  of  the 
calibration  factor  affects  R  linearly,  and  it  affects  Rk  exponentially.  This  explains  results  of 
above  Cases  1-3. 

We  now  explain  Figure  23  as  well  as  results  of  tables  9,10,11.  To  do  this,  we  use  (92)  to 
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Table  11:  Computed  values  of  Rk  for  Case  3.  Bush  is  the  calibration  target.  Bold  faced 
are  values  of  Rk  which  lead  to  values  of  Er  {target)  outside  of  tabulated  intervals. 


target 

Rk 

Er  {bckgr) 

computed  Er  {target) 

tabulated  Er  {target) 

metal  cylinder 

2.0 

|3,5] 

6,10] 

[10,30] 

metal  box 

1.16 

|3,5] 

[3.48,5.8] 

[10,30] 

wood  stake 

3.2 

1 

3.2 

|2.61 

plastic  cylinder 

0.74 

|3,5] 

[3.22,3.7] 

1.2 

70.0 

60.0 

50.0 


AGCM 

KIEM 


40.0 


(b) 


30.0 

20.0 

10.0 

0.0 


wood  stake 


0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0 


Figure  23:  Dependencies  of  R  and  Rk  from  the  calibration  factor  for  bush  and  wood  stake 
are  displayed.  In  both  a),b)  “1”  on  the  y— axis  for  AC  CM  (the  globally  convergent  method) 
corresponds  to  CF  =  10“’^  in  (90).  a)  “1”  on  the  y—axis  in  KIEM  (Krein  method)  corre¬ 
sponds  to  CFK{bush)  =  0.6-10“^  (Case  3).  b)  “F’  on  the  y—axis  in  KIEM  (Krein  method) 
corresponds  to  CFK{wood  stake)  =  10“®  (Case  1). 
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present  the  Krein  equation  (91)  as 


Z 

p{z,t) j  j'{t-T)p{z,T)dT  =  ^,  te[-z,z],  V^e[0,T/2].  (93) 

—z 


In  (93)  fit)  is  the  pre-preprocessed  data  as  in  Stage  1  and  /(+0)  =  1  as  on  Figure  23. 
Denote  /3  :=  z  ■  CFK supi^j^^  f  [t]  .  If  /5  G  (0, 1),  then  one  can  solve  integral  equation  (93) 

via  the  classical  resolvent  series.  It  is  clear  from  this  series  that  the  solution  piz,  t)  changes 
almost  with  an  exponential  speed  when  the  calibration  factor  CFK  changes. 

On  the  other  hand,  taking  into  account  the  calibration  factor  CF  for  the  globally  con¬ 
vergent  method,  the  function  q  (x,  s)  in  (18)  can  be  represented  as 


In  {CF  ■  w  {x,  s)) 


In  {CF  ■  w  {x,  s))  ^  Ws  {x,  s) 


Ws  (a;,  s)J 


Since  we  work  with  large  values  of  s  >  1,  then  q  {x,  s)  changes  logarithmically  with  the 
change  of  CF.  This  is  a  much  slower  rate  of  change  than  the  exponential  one  of  the  Krein 
equation. 


8  Global  Convergence  via  a  Carleman  Weight  Func¬ 
tion 


The  statement  of  the  CIP3.1  in  section  3.1  uses  only  the  Dirichlet  boundary  condition  (19). 
However,  using  the  procedure,  which  is  a  simple  generalization  of  the  data  complementing 
in  (76),  it  is  possible  to  prove  that  one  can  approximately  hnd  both  Dirichlet  and  Neumann 
boundary  conditions  for  the  function  q  {x,  s)  on  the  entire  boundary  dQ.  Thus,  by  (18),  we 
should  solve  the  following  problem  for  the  function  q  {x,  s) 


Aq  —  2s^Vg 


OO 

p 

\  ^  1 

1  /* 

/  Vq  {x,  t)  dr  +  2s 

<1 

Si¬ 

'S 

kJ 

s 

1  Kf 

\  -S  1 

0,  a;  G  D,  s  >  s. 


q  1^0=  {x,  s) ,  duq  |af7=  ip{x,s),s>  s. 


(94) 

(95) 


Since  the  first  work  of  Beilina  and  Klibanov  [25]  on  the  topic  of  global  convergence,  the 
approach  to  the  problem  (94),  (95),  which  has  been  pursued  so  far  and  which  is  the  focus 
of  the  effort  of  this  project,  by  the  PI  is  based  on  the  truncation  of  the  integral  in  (94)  at  a 
high  value  of  the  pseudo  frequency  s  :=s.  In  other  words,  it  has  been  assumed  so  far  that 


OO 

s 


q  {x,  t)  dr  +  V  {x,  s)  ,s  E  (s,  s) . 
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Next,  the  function  V  {x,s),  which  is  called  “the  tail  function”,  and  the  function  q  {x,  s)  were 
approximated  separately  in  an  iterative  procedure.  In  particular,  the  first  approximation 
Vi,i  (x)  for  the  function  V  {x,  s)  was  found  via  truncation  of  an  asymptotic  series  with  respect 
to  1/s,  s  >>  1  and  solution  of  the  resulting  Dirichlet  boundary  value  problem  (39),  (40)  for 
the  Laplace  equation. 

Unlike  the  above,  it  was  hgured  out  in  the  recent  work  [11]  how  to  work  out  a  globally 
convergent  numerical  method  without  any  truncation  of  the  integrals  in  (94).  Instead,  a 
signihcantly  new  idea  is  used.  Namely,  a  weighted  cost  functional  is  constructed  in  [11]. 
The  weight  is  the  Carleman  Weight  Function  (CWF),  which  is  involved  in  the  Carleman 
estimate  for  the  Laplace  operator.  Let  H  he  a  certain  Hilbert  space.  Given  a  bounded  set 
G  C  FT  of  an  arbitrary  size,  one  can  choose  a  parameter  of  CWF  in  such  a  way  that  that 
functional  becomes  strictly  convex  on  G.  The  strict  convexity,  in  turn  guarantees  the  global 
convergence  of  the  gradient  descent  method,  i.e.  this  method  can  start  from  any  point  of  G 
and  still  will  converge  to  the  exact  solution. 

Thus,  the  work  [11]  represents  a  new  globally  convergent  numerical  method  for  CIP3.1. 
Another  new  approach  to  the  issue  of  the  global  convergence  via  another  CWF  was  recently 
proposed  in  [12]. 

Globally  strictly  convex  cost  functionals  with  Carleman  Weight  Functions  in  them  were 
constructed  by  the  PI  in  his  previous  publications  of  1997  [49,  50].  However,  those  func¬ 
tionals  were  constructed  in  the  time  domain  for  CIPs  for  hyperbolic  [49]  and  parabolic  [50] 
PDFs.  Equation  (94),  which  is  in  the  Laplace  transform  domain,  was  not  considered  in  these 
references  and  the  problem  of  tails  was  not  addressed. 

8.1  Outline  of  the  new  globally  convergent  numerical  method  of 

[11] 

We  represent  the  function  q{x,  s)  as  a  series  with  respect  to  an  orthonormal  basis  of  L2{s,  oo). 
Using  this  representation,  the  integral  over  the  inhnite  interval  (s,  cxd)  in  (94)  can  be  easily 
computed.  Let  {/„  (s)})]^o  ^  -^2  (^  (d) ,  cxd)  be  an  orthonormal  basis  in  L2  (s  (d) ,  cxd)  such 
that  {fn  (s)}[]Ug  C  Li  (s  (d) ,  cx)) .  As  an  example,  one  can  consider  Laguerre  functions  [24] 

n  ^  I 

(*)  =  E  (-i)‘  ^  <“■  ’  T  = 

Next,  we  set  /n  (<s)  :=  L„(s  — s(d)),s  G  {s{d),oo).  It  can  be  verihed  that  q{x,s)  G 
L2  (s  (d) ,  cxd)  ,  Vx  G  C.  Hence,  one  can  represent  the  function  q  {x,  s)  as 

OO  A^— 1 

Q  (a^,  s)  =  ^qn  (x)  fn  (s)  ^  ^  g„  (a;)  fn{s),  s>s  (d) ,  (96) 

n=0  n=0 

where  iV  is  a  sufficiently  large  integer  which  should  be  chosen  in  numerical  experiments. 
Consider  the  vector  of  coefficients  in  the  truncated  series  (96)  Q  {x)  =  (go,  •••,  Qn-i)  {x)  G 
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Substituting  the  truncated  series  (96)  in  (94),  we  obtain 


N-l 


N~l  N-l 


^Qn{x)Us)  m{,x)Vqn{x)f^{s) 

n=0  m=0  n=0 


oo 


s 


fn{r)dT 


N-l  N-l 

+  2s  ZE  qmix)'\/ Qnix) 
m=0  n=0 


fm{r)dT  /  fn{r)dT 


0. 


(97) 


To  be  precise,  one  should  have  instead  of  “=”  in  (97)  due  to  the  truncation  (13). 
Multiplying  both  sides  of  (97)  by  fk{s),  integrating  over  {s{d) ,  cxd)  and  keeping  in  mind  the 
fact  that  {/n('S)}))^o  orthonormal  basis  in  ^2(5,  cx)),  we  obtain  the  following  system  of 
coupled  nonlinear  elliptic  equations: 

N-l  N-l 

Aqk{x)  +  FkmnVqm{x)Vqn{x)  =  0,  k  =  0,...,N -1,  X  eQ,  (98) 

m=0  n=0 

where  the  numbers  F^mn,  k,m,n  E  {0, . . . ,  iV  —  1},  are  given  by 

OO  OD  OO  OO  OO 

Fkmn=  J  2sfkis){  j  fm{T)dT  J  /„(r)dr)ds-  J  2s'^fkis)fmis){  j  fn{,r)dT)ds. 

s{d)  S  s  s(d)  S 


The  boundary  conditions  for  qn  are  obtained  by  substituting  again  the  truncated  series 
(96)  into  (95).  For  the  convenience  of  the  following  analysis,  we  rewrite  system  (98)  together 
with  the  boundary  conditions  as  the  following  boundary  value  problem  with  over-determined 
boundary  conditions.  Note  that  we  have  both  Dirichlet  and  Neumann  boundary  conditions 


AQ  +  F(VQ)  =  0, 

Q  |so=  (x) ,  dnQ  |ao=  T  (x) , 


(99) 

(100) 


where  the  boundary  vector  functions  <h  (x) ,  T  (x)  E  are  computed  from  the  functions 
if  (x,  s) ,  ^  (x,  s)  and  F  :  ^  R^,  F  =  (Fq,  ...,  F^_i)  e  (R^^)  with 


Af-l  N-l 


Fk{VQ)  =  EE  Fkmn^ qm (x)^ qn (x^ ,  k  0,  .  .  .  ,  Af  1. 

m=0 n=0 

If  we  can  hnd  an  approximate  solution  of  the  problem  (99)-(100),  then  we  can  hnd  an 
approximation  for  the  function  q  via  the  truncated  series  (96).  Therefore,  we  focus  below  on 
the  method  of  approximating  the  vector  function  Q  [x) . 

Let  F  (x)  =  (Fq,  ...,  Fat-i)  (x)  ,x  E  kl,  be  a  vector  function  and  F  be  a  Hilbert  space. 
Below  any  statement  that  F  E  H  means  that  every  component  of  the  vector  F  belongs  to 


H.  The  norm 


means 


H 


F 


H 


=  ( 


N-l 

E 

n=0 


H/2 


H 
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Our  ultimate  goal  is  to  apply  this  method  to  the  inversion  of  above  experimental  data 
of  [2,  3,  4,  5].  Thus,  just  as  in  these  references,  below  O  is  chosen  to  be  a  rectangular 
parallelepiped.  Without  loss  of  generality,  it  is  convenient  to  assume  that 

O  =  {x  =  {xi,X2,X3)  :  {xi,X2)  e  {-A,  A)  ,X3  e  (0, 1/2)} , 

where  A  >  0  is  a  number.  Thus,  dQ  =  riUr2Ur3,  where 

Ti  =  {x  e  dn\xs  =  0}  ,  Ts  =  {x  e  dn\xs  =  1/2}  ,  Ts  =  0\  {T,  U  Ta) . 

As  in  [2,  3,  4,  5],  Ti  is  considered  as  the  backscattering  side,  where  the  data  are  measured.  Al¬ 
though  measurements  were  not  performed  on  r2Ur3,  it  was  demonstrated  in  these  references 
that  assigning 

\r^ur3-=  Woix,s)  |r2ur3  (101) 

does  not  affect  the  accuracy  of  the  reconstruction  via  the  technique  of  [1].  Thus,  we  now 
relax  conditions  (100),  assuming  that  the  normal  derivative  is  given  only  on  Ti,  the  Dirichlet 
condition  is  given  on  Ti  U  Ts  and  no  boundary  condition  is  given  on  r2. 


Q  IriUFs—  ^  )  dnQ  Iti—  d'  (x) . 


(102) 


Let  us  introduce  a  CWF  for  the  Laplace  operator  which  is  suitable  for  this  domain  and 
for  boundary  conditions  (102).  Let  a,  ^  G  (0,1/2)  be  two  arbitrary  numbers.  Let  A,  z/  >  1 
be  two  large  parameters  which  we  will  choose  later.  Then  the  CWF  has  the  form 

(xg)  =  e^C3+«)-"e-A(a+0-"_  (103) 

Hence, 

lim  ^x,u  (1/2)  =  0. 

A— >-oo 

Lemma  9.1  establishes  a  Carleman  estimate  for  the  operator  A  in  the  domain  with  the 
weight  function  (103). 

Lemma  9.1.  There  exist  sufficiently  large  numbers  Aq  =  Aq  (fl)  >  1,  z/q  =  (H,  0  >  1 

depending  only  on  the  listed  parameters  such  that  for  an  arbitrary  function  u  G  (12) 
satisfying  u  |r3=  0  the  following  Carleman  estimate  holds  for  all  A  >  Aq  and  with  a  constant 
C  =  C  (12)  >  0  depending  only  on  the  domain  12 

J  +  C  (^\\u  Irill^pri)  +  IrillL(ri)) 

n 

>C  j  (AIVuIVaV)  j  {\V  u\^  +  ffi)  dx2dx:,. 


Below  C  =  C  (12)  >  0  denotes  different  constants  depending  only  on  the  domain  12.  Let 
i?  >  0  be  an  arbitrary  number.  Dehne  the  set  G  of  vector  functions  Q  as 


Q  =  (go, gv-i)^  e  (12)  :  ||g||^3(^)  < 
Q  |riur3=  d>  (x) ,  d„Q  |ri=  'L  (x) . 


Then  G  is  an  open  set  in  (hi) .  Also,  Embedding  Theorem  implies  that 

GcC‘ (SI),  ||(J||p,p)<Cfl,VOeG. 

Let  uq  =  (n,  a,  .^)  >  1  be  the  number  in  Lemma  9.1.  Denote  =  Dfllxs  <  a}  .  We  seek  the 
solution  Q  of  the  problem  (99),  (102)  on  the  set  G  via  minimizing  the  following  Tikhonov-like 
cost  functional  with  the  CWF  and  with  the  regularization  parameter  a  G  (0, 1) 

Jx.a  (Q)  =  \j  |AO  +  F  (VQ)|"  ^fl„dx  +  I  ,  Q  €  G  (fl,  4, 4-) . 

o 

Theorem  9.2.  There  exists  a  sujficiently  large  number  Ai  =  Ai  {Q,G,F)  >  1  depending 
only  on  il,G,F  such  that  if  \  >  \i  and  a  E  functional  J\a{Q) 

is  strictly  convex  on  the  set  G,  i.e.,  there  exists  a  constant  Gi  =  Gi  {Q,G,F)  >  0  depending 
only  on  D,  G,  F  such  that  for  all  Qi,Q2  ^  G 

j\,a  {Q2)  —  Jx,a  (Ql)  —  J'x^a  (Qi)  (Qs  “  Ql) 

>  Cl  \\Q2  -  Qi\\%i(q^)  +  —  IIQ2  -  Qi|Ih3(0)  > 

where  (Qi)  is  the  Frechet  derivative  of  the  functional  Jx^a  at  the  point  Qi. 

8.2  Global  convergence  of  the  gradient  descent  method 

It  is  well-known  that  the  gradient  descent  method  is  globally  convergent  for  functionals 
which  are  strictly  convex  on  the  entire  space.  However,  the  functional  Ja,q  is  strictly  convex 
only  on  the  bounded  set  G(i?,  <h,T).  Therefore,  we  need  to  prove  the  global  convergence 
of  this  method  on  this  set.  Suppose  that  a  minimizer  Qmin  of  Ja,«  exists  on  G(i?,  <h,  T). 
In  the  regularization  theory  Qmin  is  called  regularized  solution  of  the  problem  (99),  (102) 
[1].  Theorem  9.2  guarantees  that  such  a  minimizer  is  unique.  First,  we  estimate  in  this 
section  the  distance  between  regularized  and  exact  solutions,  depending  on  the  level  of  error 
in  the  data.  Next,  we  establish  that  Theorem  9.2  implies  that  the  gradient  descent  method 
of  the  minimization  of  the  functional  Jx^a  converges  to  Qmin  if  starting  at  any  point  of  the 
set  G  (i?/4,  $,  T),  i.e.,  it  converges  globally.  In  addition,  we  estimate  the  distance  between 
points  of  the  minimizing  sequence  of  the  gradient  descent  method  and  the  exact  solution  of 
the  problem.  In  principle,  global  convergence  of  other  gradient  methods  for  the  functional 
Jx,a  can  also  be  proved.  However,  we  are  not  doing  this  for  brevity. 

Following  one  of  concepts  of  Tikhonov  for  ill-posed  problems  (see,  e.g.,  section  1.4  in  [1]), 
we  assume  that  there  exist  noiseless  boundary  data  {x)  and  T*  [x)  which  correspond  to 
the  exact  solution  Q*  of  the  problem  (99),  (102).  Also,  we  assume  that  functions  (x)  and 
4/  [x)  at  the  part  Ti  of  the  boundary  contain  an  error  of  the  level  <5, 
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On  the  other  hand,  we  do  not  assume  any  error  in  the  function  $  at  r2  U  Fs,  see  a  heuristic 
condition  (101),  which  was  justihed  numerically  in  [2,  3,  4,  5].  Theorem  9.3  estimates 
the  distance  between  Qmin  and  Q*  in  the  norm  of  (0^) ,  which  might  be  sufficient  for 
computations.  Note  that  while  in  Theorem  9.2  we  have  compared  functions  Qi  and  Q2 
satisfying  the  same  boundary  conditions,  functions  Qmin  and  Q*  in  Theorem  9.3  satisfy 
different  boundary  conditions,  because  of  the  error  in  the  data. 

Theorem  9.3.  Assume  that  conditions  of  Theorem  9.2  hold  and  \  >  Xi-  In  addition, 
assume  that  conditions  (IO4)  are  satisfied  and  also  <h  |r3=  Its-  Suppose  that  there  exists  an 
exact  solution  Q*  of  the  problem  (99),  (102)  and  Q*  G  G  (i?,  $*,  T*) .  In  addition,  assume 
that  there  exists  a  minimizer  Qmin  G  of  the  functional  J\^a-  Lot  the  number 

^0  =  ^0  G,  F,  a,  f)  G  (0, 1)  be  so  small  that 


5, 


-^-^0/2 


>  Ai. 


Let  5  G  (0,  (5o) .  Choose  the  regularization  parameter  a  as  a  =  a  (5)  =  where 

r  r,  r 


2  (a  +  (r  L‘  (“  +  0"“ 

Then  a  G  {t\,uo  (1/2)  >  l)  (os  in  Theorem  9.2)  and 


0,- 

2 


IIQ*  —  Qmin||Hl(0„)  ^  Ci6^,  WS  G  (0,  5o)  . 

We  now  formulate  the  gradient  descent  method  with  the  constant  step  size  fl  for  the 
problem  of  the  minimization  of  the  functional  J\^a-  For  brevity  we  do  not  indicate  the 
dependence  of  functions  Qn  on  parameters  A,  a,  (3.  Let  Qi  E  G  (i?/4,  <h,  T)  be  an  arbitrary 
point  of  the  set  G  (i?/4,  <h,  T).  Consider  the  sequence  of  the  gradient  descent 

method, 

Qn+l  =Qn-  /3Jx,a  (Qn)  ,  n  =  1,  2,  ...  (105) 

Theorem  9.4  establishes  the  global  convergence  of  the  sequence  (105)  on  the  set  G  {R,  <h,  T) . 
In  addition,  it  estimates  the  distances  between  points  of  this  sequence  and  the  exact  solution, 
depending  on  the  level  <5  of  the  noise  in  the  data  in  (104). 

Theorem  9.4.  Choose  parameters  Ai,z/o,q;  as  in  Theorem  9.2  and  let  A  >  Ai.  As¬ 
sume  that  the  functional  Jx^a  achieves  its  minimal  value  on  the  set  G(i?,  <h,4/)  at  a  point 
Qmin  G  G  (i?/4,  $,  T).  Consider  the  sequence  (105),  in  which  the  starting  point  Qi  is  an 
arbitrary  point  of  the  set  G  (i?/4,  $,  T).  Then  there  exists  a  sufficiently  small  number 
(3  =  /9  (A,  a,  G  (/?,  <h,  4/))  G  (0,1)  and  a  number  9  =  9  {(3)  G  (0,1),  both  dependent  only 
on  the  listed  parameters,  such  that  the  sequence  (105)  converges  to  the  point  Qmin, 


||f?n+l  Qmin  1 1  ^  IIQl  Qmin  || -,01  1,  2,  ... 

In  addition,  assume  that  there  exists  an  exact  solution  Q*  E  G  {R/5,  <h*,  T*)  of  the  problem 
(14),  (45)  and  that  all  conditions  of  Theorem  9.3  hold.  Then  for  all  S  E  (0,5o)  the  following 
estimate  holds 


IIQn+l  —  Q*ll//1(0„)  <  IIQl  ~  Qmin||//3(S2)  +  Gi^'*',  U  — 


=  1,2, 


70 


8.3  Numerical  examples 

In  this  subsection  we  present  a  limited  testing  of  the  above  algorithm  for  some  numerically 
simulated  data.  We  also  compare  its  performance  with  the  above  locally  convergent  method 
alone  using  the  coefficient  of  the  homogeneous  medium  as  the  first  guess.  Numerical  results 
for  experimental  data  in  both  1-d  and  3-d  cases  are  under  consideration  and  will  be  reported 
in  future  work. 
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Figure  24:  Reconstruction  of  the  coefficient  of  Example  1.  (a)  Result  of  Steps  1  and  2,  (b) 
Result  of  Step  3,  (c)  Result  of  Step  2  starting  from  the  homogeneous  medium  as  the  first 
guess,  (d)  Result  of  Step  3  applied  to  the  result  of  (c). 


Since  our  target  application  is  in  imaging  of  an  abnormal  object  placed  in  a  homoge¬ 
neous  medium,  we  mainly  test  the  proposed  algorithm  with  a  discontinuous  coefficient.  The 
locations  of  the  discontinuities  represent  the  location  of  the  target.  It  is  hard  to  obtain 
accurate  reconstructions  at  locations  far  away  from  the  measurement  point.  However,  this 
is  not  a  serious  restriction  from  the  practical  standpoint.  Indeed,  our  experience  of  working 
with  3-d  time  resolved  backscattering  experimental  data  [2,  3,  4,  5]  tells  us  that,  using  the 
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so-called  data  propagation  procedure  in  data  pre-processing,  we  can  approximate  quite  well 
both  the  distance  from  the  measurement  point  to  the  target  and  the  Dirichlet  and  Neumann 
backscattering  data  near  the  target.  Thus,  we  assume  below  that  the  target  is  close  to  the 
measurement  point. 

In  the  following  examples,  the  parameters  were  chosen  as  follows:  The  pseudo  frequencies 
were  s  G  [4, 15]  with  the  integration  step  size  in  (97)  As  =  0.05.  The  number  of  Laguerre’s 
functions  was  N  =  11.  The  coefficients  A  in  the  CWF  was  A  =  3.  The  regularization 
parameter  a  =  0.001  and  the  truncation  value  e  =  0.2. 
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Figure  25:  Reconstruction  of  the  coefficient  of  Example  2.  (a)  Result  of  Steps  1  and  2,  (b) 
Result  of  Step  3. 


Example  1.  Consider  a  piecewise  constant  coefficient  given  by  c{x)  =  1  -|-  3x[0.03,  0.1], 
where  y  is  the  characteristic  function.  Figure  24  (a)  compares  the  computed  coefficient 
values  of  Steps  1  and  2  with  the  exact  one  and  Figure  24  (b)  depicts  the  result  of  Step  3. 
To  compare  with  the  performance  of  the  above  locally  convergent  method  alone,  we  plot  in 
Figure  24  (c),  (d)  the  results  of  Steps  2  and  3,  respectively,  starting  from  the  homogeneous 
medium  as  the  initial  guess.  We  can  see  that  our  hybrid  algorithm  provided  accurate  results, 
whereas  the  locally  convergent  method  alone  failed. 

Example  2.  In  this  example,  we  consider  another  piecewise  constant  coefficient  with  a 
larger  jump  c{x)  =  1  -|-  14y;[0.03,  0.1].  The  results  of  Steps  1-3  are  shown  in  Figure  25.  Even 
though  the  jump  of  the  coefficient  is  high  in  this  case,  we  still  can  see  the  good  accuracy  of 
the  reconstruction  using  our  hybrid  algorithm. 

Example  3.  Consider  the  exact  coefficient  c{x)  =  1  —  0.5x[0.03,  0.15].  This  coefficient 
mimics  the  case  the  case  when  the  dielectric  constant  of  an  explosive  is  less  than  the  one  of 
a  homogeneous  background.  Figure  26  shows  the  reconstruction  results  of  Steps  1  -  3  of  the 
algorithm.  Again,  we  obtained  an  accurate  reconstruction. 

Example  4.  Finally,  we  consider  a  continuous  coefficient  given  by  c{x)  =  /(0-04)  _ 

The  results  are  shown  in  Figure  27.  Comparing  Figure  27  (a)  with  Figure  27  (c),  one  can  see 
that  the  combination  of  Step  1  and  Step  2  provided  much  better  result  than  Step  2  starting 
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X  X 

(a)  (b) 

Figure  26:  Reconstruction  of  the  coefficient  of  Example  3.  (a)  Result  of  Steps  1  and  2,  (b) 
Result  of  Step  3. 

from  the  homogeneous  medium  as  the  hrst  guess.  However,  results  of  Step  3  of  both  cases 
are  accurate. 

From  Figure  24  (d)  and  Figure  27  (d)  we  see  that  the  above  locally  convergent  method, 
taking  alone,  is  unstable  in  the  sense  that,  depending  on  the  type  of  the  target,  it  provides 
either  bad  or  good  quality  images.  Meanwhile,  the  proposed  hybrid  algorithm  provides 
accurate  results  in  all  four  examples. 
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(c) 


(d) 


Figure  27:  Reconstruction  of  the  coefficient  of  Example  4.  (a)  Result  of  Steps  1  and  2,  (b) 
Result  of  Step  3,  (c)  Result  of  Step  2  starting  from  the  homogeneous  medium  as  the  hrst 
guess,  (d)  Result  of  Step  3  applied  to  the  result  of  (c). 
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9  The  First  Solution  of  a  Long  Standing  Problem  [8, 
9,  10] 

In  the  paper  [8]  the  PI  has  solved,  for  the  first  time,  a  long  standing  problem  posed  in  1977 
by  French  mathematicians  Chadan  and  Sabatier  in  their  famous  book  [38].  However,  only 
after  the  paper  [8]  was  published,  the  PI  has  realized  that  the  result  of  this  paper  indeed 
represents  the  first  solution  of  that  problem.  Thus,  the  PI  has  published  a  short  paper  [9], 
where  he  has  pointed  out  to  this  fact.  A  generalization  of  the  result  of  [8]  to  the  case  of  the 
acoustics  equation  in  was  published  in  the  paper  of  the  PI  [10]. 

When  considering  Coefficient  Inverse  Problems  in  the  frequency  domain,  it  was  assumed 
in  all  previous  publications  that  both  the  modulus  and  the  phase  of  the  complex  valued 
function  representing  the  wave  held  is  known  on  a  certain  set,  see,  e.g.  [60,  61]  for  global 
uniqueness  results  and  reconstruction  methods.  However,  in  quantum  inverse  scattering, 
which  is  the  main  interest  of  the  book  [38],  only  the  differential  cross-section  is  measured. 
This  means  that  only  the  modulus  of  the  scattering  wave  is  measured  in  quantum  inverse 
scattering.  On  the  other  hand,  the  entire  theory  of  the  quantum  inverse  scattering  uses  the 
assumption  that  both  the  modulus  and  the  phase  of  the  scattering  complex  valued  wave  held 
are  measured.  This  caused  the  authors  of  the  book  [38]  to  pose  in  1977  in  chapter  10  the 
question  about  the  uniqueness  of  the  inverse  scattering  problem  in  the  case  when  only  the 
modulus  of  the  scattering  complex  valued  wave  held  is  measured  outside  of  the  scatterer.  As 
it  was  mentioned  above,  this  question  was,  for  the  hrst  time,  addressed  by  the  PI  in  2014  in 
[8].  Unlike  the  current  3-d  case,  in  the  past,  uniqueness  of  the  phaseless  Coefficient  Inverse 
Problem  was  proven  only  for  the  case  of  the  1-d  Schrodinger  equation  by  Klibanov  and  Sacks 
in  1992  [53]. 

Inverse  problems  without  the  phase  information  are  well  known  in  optics,  since  it  is  often 
impossible  to  measure  the  phase  of  the  optical  signal,  unlike  its  amplitude.  In  optics,  such  a 
problem  is  usually  formulated  as  the  problem  about  the  recovery  of  a  compactly  supported 
complex  valued  function  from  the  modulus  of  its  Fourier  transform.  The  latter  is  called  the 
“phase  retrieval  problem”  [43,  44].  Uniqueness  theorems  for  the  phase  retrieval  problem  can 
be  found  in  [52,  54,  55]. 

9.1  The  first  solution  of  the  Chadan-Sabatier  problem 

Below  are  Holder  spaces,  where  s  >  0  is  an  integer  and  a  G  (0, 1) .  Let  H,  G  C  R^  be 
two  bounded  domains,  Q  G  G.  For  an  arbitrary  point  ?/  G  R^  and  for  an  arbitrary  number 
u  G  (0, 1)  denote  {y)  =  {x  :  \x  —  y\  <  u}  and  {y)  =  R^\Pt^  (y) .  For  any  two  sets 
M,  A^  C  R^  let  dist  (M,  N)  be  the  Hausdorff  distance  between  them.  Let  Gi  C  R^  be  a 
convex  bounded  domain  with  its  boundary  S  G  G^.  Let  e  G  (0, 1)  be  a  number.  We  assume 
that  H  C  Gi  C  G,  dist  (S',  dG)  >  2e  and  dist  (S',  (912)  >  2e.  Hence, 

dist  {dB^  (y) ,  d^l)  >  e,Vy  E  S', 

dist  {dBs  (y) ,  dG)  >  e,'iy  E  S'. 
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(106) 

(107) 


Below  either  m  =  2  or  m  =  4,  and  we  will  specify  this  later.  We  impose  the  following 
conditions  on  the  potential  q  {x) 


q  (x)  e  C™  (M^)  ,  g  (x)  =  0  for  a;  e  (108) 

q{x)  >  0.  (109) 

Let  xq  =  (xo,!,  a;o,2, 2:0,3)  be  the  sonrce  position.  Consider  the  following  problem 

^xU  +  k‘^u  —  q  {x)  u  =  —5  {x  —  xq)  ,  x  G  R^,  (110) 


- ^dx  u  (x,  Xq,  k)  —  iku  (x,  Xq,  k)  =  O  (1)  ,  Ixl  — )■  CXD. 

“  \x-Xo\  " 

J=1 


(111) 


Here  the  radiation  condition  (111)  is  valid  for  every  hxed  sonrce  position  xq.  To  establish 
existence  and  nniqueness  of  the  solution  of  the  problem  (llO)-(lll),  we  refer  to  Theorem  6  of 
Chapter  9  of  the  book  [68]  as  well  as  to  Theorem  3.3  of  the  paper  [67].  As  to  the  smoothness 
of  the  solution  of  the  problem  (llO)-(lll),  we  refer  to  Theorem  6.17  of  the  book  [45].  Thus, 
combining  these  results,  we  obtain  that  for  each  pair  {k,Xo)  G  R  x  R^  there  exists  unique 
solution  u  {x,Xo,  k)  of  the  problem  (110),  (111)  such  that 


U  {x,  Xq,  k)  =  Mo  {.X,  Xq,  k)  +  Ug  (x,  Xq,  k)  ,  (112) 

^  exp(^Xa:-:ro|)^  ^  ^ ^  ^  _ 

47r  \x  —  xo\ 

In  (112),  (113)  uo{x,xo,k)  is  the  incident  spherical  wave  and  Us{x,xo,k)  is  the  scattered 
wave. 

Phaseless  Inverse  Problem  1  (PIPl).  Let  m  =  2  in  (108).  Suppose  that  the  function 
q{x)  satisfying  (108),  (109)  is  unknown  for  x  E  kl  and  known  for  x  G  R^\f2.  Also,  assume 
that  the  following  function  fi  {x,Xo,  k)  is  known 


fi  {x,Xo,k)  =  \u{x,Xo,k)  \  ,Va:o  E  S,Wx  e  {xq),x  ^  XQ,\/k  E  {a,b) , 


where  (a,  &)  C  R  is  an  arbitrary  interval.  Determine  the  function  q  [x)  for  x  E  kl. 

We  now  outline  the  main  difficulty,  which  did  not  allow  one  to  prove  uniqueness  results  for 
phaseless  3-d  inverse  scattering  problems  so  far.  As  an  example,  consider  IPl.  Analogous 
difficulties  take  place  all  other  phaseless  Coefficient  Inverse  Problems  formulated  in  this 
report.  In  PIPl  one  should  work  with  a  complex  valued  function  r  (A;) ,  A;  G  R  such  that  its 
modulus  |r  [k)  \  is  known  for  all  k  E  (a,  h).  The  function  r  (A;)  admits  the  analytic  continuation 
from  the  real  line  R  in  the  half-plane  {A;  G  C  :  Im  A;  >  —7}  for  a  certain  number  7  >  0.  Since 
|r  {k)^  =  r  {k)f  [k) ,  then  the  function  |r  {k)\^  is  analytic  for  A;  G  R  as  the  function  of  the 
real  variable  k.  Here  r  {k)  is  the  complex  conjugate  of  r  (A;) .  Hence,  the  modulus  |r  (A;)|  is 
known  for  all  A;  G  R.  Denote  C+  =  {A;  G  C  :  ImA;  >  0} .  Proposition  4.2  of  [54]  implies  that 
if  r  (A;)  would  not  have  zeros  in  C+,  then  this  function  would  be  uniquely  reconstructed  for 
all  A;  G  R  from  the  values  of  |r(A;)|  for  A:  G  R.  However,  the  main  dijficulty  is  to  properly 
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account  for  zeros  of  r  {k)  in  the  upper  half-plane  C+\M.  Indeed,  let  zi,  ■■■,Zn  G  C+\M  be 
some  of  such  zeros  of  r  [k) .  Consider  the  function  r  (/c)  dehned  as 

Hence,  |r  (/c)|  =  \r  {k)  \ ,  V/c  G  M.  Furthermore,  the  function  r  (/c)  is  analytic  in  the  half-plane 

{k  E  C  :lmk  >  —7}  .  Therefore,  in  order  to  prove  uniqueness,  one  needs  to  hgure  out 

how  to  combine  the  knowledge  of  |r  {k)\  for  A;  G  M  with  a  linkage  between  the  function  r  [k) 

and  the  originating  differential  operator. 

This  difficulty  was  handled  in  [53]  in  the  1-d  case,  using  the  fact  that  the  function  r  [k) 
depends  only  on  one  variable  k  in  this  case.  Unlike  [53],  the  function  r  (x,  Xq,  k)  =  u  {x,  Xq,  k) 
depends  on  x,Xq,  k  in  the  3-d  case.  Hence,  the  above  zeros  depend  now  on  both  x  and  Xq, 
i.e.  Zj  =  Zj  (x,  xo) .  Thus,  compared  with  the  1-d  problem,  the  main  difficulty  of  the  3-d 
case  is  that  it  is  necessary  to  hgure  out  how  to  take  into  account  the  dependence  of  zeros 
Zj  {x,Xo)  from  x  and  Xq.  To  do  this,  we  essentially  use  here  properties  of  the  solution  of  the 
Cauchy  problem  for  an  associated  hyperbolic  PDF 

Vtt  =  An  —  q  (x)  n,  {x,  t)  G  x  (0,  cxd)  ,  (114) 

V  {x,  0)  =  0,  Vt  {x,  0)  =  (5  (x  —  xo) .  (115) 

Indeed,  it  follows  from  Lemma  6  of  Chapter  10  of  the  book  [68]  as  well  as  to  Remark  3  after 
that  lemma  that  the  function  u  is  the  Fourier  transform  of  the  function  v, 

00 

u{x,Xo,k)  =  j  V  {x,  Xo,t)  Vx,  Xq  G  a;  7^  Xq,  V/c  G  M.  (116) 

0 

Theorem  9.1.  Consider  PIPl.  Let  two  potentials  qi  {x)  and  q2  {x)  satisfying  condi¬ 
tions  (108),  (109)  be  such  that  qi{x)  =  q2{x)  =  q{x)  for  x  G  M^\r2.  Let  ui{x,xo,k) 
and  U2  {x,Xf),  k)  he  corresponding  solutions  of  the  problem  (llO)-(lll)  satisfying  conditions 
(112),  (113).  Assume  that 

\ui  {x,  Xo,k)  \  =  \u2  {x,  Xo,k)  \ ,  Vxo  G  S',  Vx  G  (xq)  ,  a;  7^  Xq,  V/c  G  (a,  5) .  (117) 

Then  qi  {x)  =  q2  {x) . 

In  PIPl  the  modulus  of  the  total  wave  held  n  =  wq  +  is  known  on  a  certain  set,  see 
(112),  (113)  for  uo  and  Us-  We  now  consider  the  case  when  only  the  modulus  of  the  scattered 
wave  is  known. 

Phaseless  Inverse  Problem  2  (PIP2).  Let  m  =  2  in  (108).  Suppose  that  the  function 
q{x)  satisfying  (108),  (109)  is  unknown  for  x  E  kl  and  known  for  x  E  M^\r2.  Also,  assume 
that  the  following  function  /2  (a;,  xq,  k)  is  known 

/2  (a;,  xo,  k)  =  \us  {x,  XQ,k)\,  Vxq  E  S,Vx  E  Be  (xq)  ,  x  7^  xq,  V/c  G  (a,  b) . 
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Determine  the  function  q  [x)  for  x  E  Q. 

Theorem  9.2.  Consider  IP2.  Assume  that  all  conditions  of  Theorem  5.1  hold,  except 
that  (in)  is  replaced  with 

\us,i  {x,  Xq,  k)  I  =  \us,2  {x,  Xq,  k)  I ,  Vxo  e  S,  Vx  e  (xq)  ,x  ^Xo,yk  e  (a,  b) ,  (118) 

where  Ugj  =  Uj  —  uo,j  =  1,2.  In  addition,  assume  that  q{x)  7^  0,  Vx  G  S.  Then  qi  {x)  = 

Q2  {x)  . 

Theorems  5.1  and  5.2  are  formulated  only  for  the  over-determined  data.  Indeed,  in 
both  PIPl  and  PIP2  the  number  of  free  variables  in  the  data  exceeds  the  number  of  free 
variables  in  the  unknown  coefficient.  The  reason  of  this  is  that  even  if  the  phase  is  known, 
still  all  current  uniqueness  results  for  3-d  inverse  scattering  problems  in  the  case  when  the 
(5— function  is  the  source  function  are  valid  only  if  the  data  are  over-determined  ones,  see, 
e.g.  [60,  61]  for  the  frequency  domain  and  §1  of  chapter  7  of  [56]  for  an  inverse  scattering 
problem  in  the  time  domain. 

Suppose  now  that  the  function  5  {x  —  Xq)  in  (110)  is  replaced  with  such  a  function  p  (x) 
that  p{x)  ^  0  in  D.  And  consider  the  inverse  problem  of  the  reconstruction  of  the  potential 
q  {x)  from  values  of  the  function  u  {x,  k)  for  all  a;  G  S',  A;  G  M.  Then  uniqueness  theorem  for 
this  problem  can  be  proved  for  the  non-overdetermined  case.  This  proof  can  be  handled  by 
the  method,  which  was  introduced  in  the  originating  paper  [37].  Also,  see,  e.g.  [13]  about 
this  method.  This  technique  is  based  on  Carleman  estimates. 

Consider  the  function  y  (x)  G  C°°  (M^)  such  that  y  (x)  =  1  in  Gi  and  y  (x)  =  0  for  x  ^  G. 
Let  Xq  G  S.  For  a  number  u  >  0  consider  the  function  {x  —  Xq)  , 

w  ^  ry  x{x) 

da  [X  —  Xn)  =  C - O  exp 

(2V^)’ 

where  the  number  C*  >  0  is  such  that 

G 

The  function  {x  —  xq)  approximates  the  function  6  {x  —  xo)  in  the  distribution  sense  for 
sufficiently  small  values  of  a.  The  function  {x  —  xo)  is  acceptable  in  Physics  as  a  proper 
replacement  of  5  (a;  —  a;o),  since  there  is  no  “true”  delta-function  in  the  reality.  On  the  other 
hand,  the  above  mentioned  method  of  [37]  is  applicable  to  the  case  when  6  (x  —  xq)  is  replaced 
with  da  (x).  The  function  da  (x  —  Xq)  is  acceptable  in  Physics  as  a  proper  replacement  of 
d  (x  —  Xq),  since  there  is  no  “true”  delta-function  in  the  reality.  On  the  other  hand,  the 
above  mentioned  method  of  [37]  is  applicable  to  the  case  when  d  (x  —  Xq)  is  replaced  with 
da  (x).  Therefore,  it  seems  to  be  worthy  from  the  Physics  standpoint  to  consider  Inverse 
Problems  3,4  below. 

Let  in  (108)  m  =  4.  To  apply  results,  which  follow  from  the  method  of  [37],  consider  the 
function  g  (x)  such  that 

gEG^  (M^)  ,  c/  (a;)  =  0  in  (119) 
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g{x)^Q  in  Gi. 


Consider  the  following  problem 

/S.W  +  k'^w  —  q  {x)  w  =  —g  {x) ,  X  E 


(120) 

(121) 


o(l),|a;|  — )-cxD.,|a;|  -E  oo. 


(122) 


The  same  results  of  [45,  67,  68]  as  ones  cited  above  in  this  section  guarantee  that  for  each 
k  E  M.  there  exists  unique  solution  v{x,k)  E  (M^) ,  Va  E  (0,1)  of  the  problem  (121), 
(122). 

Phaseless  Inverse  Problem  3  (PIPS).  Let  m  =  4  in  (108).  Suppose  that  the  function 
q{x)  satisfying  conditions  (108),  (109)  is  unknown  for  x  E  kl  and  known  for  x  G  M^\r2. 
Assume  that  the  following  function  /s  [x,  k)  is  known 


/s  {x,  k)  =  \w  {x,  k)\  ,\/x  E  S,Wk  E  (a,  b)  . 


(123) 


Determine  the  function  q  [x)  for  x  E  kl. 

Theorem  9.3.  Consider  PIPS.  Let  the  function  g{x)  satisfies  conditions  (119),  (120). 
Consider  two  functions  qi  {x)  ,q2  (x)  satisfying  conditions  (108),  (109)  and  such  that  qi  {x)  = 
Q2  (x)  =  q  (x)  for  x  E  M^\D.  For  j  =  1,2  let  Wj  (x,  k)  E  (M^^  ^/^g  solution  of  the 

problem  (121)-(122)  with  q{x)  =  qj  {x).  Assume  that 

|tci  {x,  k) I  =  \w2  {x,  k)\  ,'ix  E  S,'ik  E  (a,  b) .  (124) 


Then  qi  {x)  =  q2  (x) . 

We  now  pose  an  analog  of  PIP2.  Let  wq  (x,  k)  be  the  solution  of  the  problem  (121)-(122) 
for  the  case  q  (x)  =0, 

f  exp  (ik  jx  -  fj) 

^)  =  y  A'K\x-f\  ^ 

G 

Hence,  one  can  interpret  the  function  wq  {x,  k)  as  the  solution  of  the  problem  (121)-(122)  for 
case  of  the  background  medium. 

Phaseless  Inverse  Problem  4  (PIP4).  Let  m  =  4  in  (108).  Suppose  that  the 
function  q{x)  satisfying  (108),  (109)  is  unknown  for  x  E  Ll  and  known  for  x  E  Let 

Wg  {x,  k)  =  w  {x,  k)  —  Wq  {x,  k) .  Assume  that  the  following  function  /4  {x,  k)  is  known 


/4  {x,  k)  =  \ws  {x,  k)\  ,Wx  E  S,Wk  E  (a,  b) . 


Determine  the  function  q  [x)  for  x  E  Ll. 

Theorem  9.4.  Consider  PIP4.  Let  all  conditions  of  Theorem  9.3  hold,  except  that  (124) 
is  replaced  with 

{x,  k)  I  =  \ws,2  {x,  /c)  I ,  Vx  G  S',  V/c  G  (a,  b) ,  (125) 

where  Wsg{x,k)  =  Wj{x,k)  —wo{x,k),j  =  1,2.  Assume  that  q{x)  7^  0,Va;  G  S.  Then 
qi  (x)  =  q2  {x) . 
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9.2  Uniqueness  of  some  phaseless  coefficient  inverse  problems  for 
the  acoustic  equation  in  3D 

We  now  formulate  uniqueness  theorems  of  the  paper  [10]  for  some  phaseless  inverse  problems 
for  the  3D  acoustics  equation.  The  main  difficulty  of  proofs  of  these  theorems  is  the  same  as 
the  one  indicated  in  section  5.3.2  for  the  case  of  the  Schrodinger  equation.  We  keep  notations 
of  the  previous  section  for  domains  and  surfaces. 

Let  c  {x)  be  the  variable  sound  speed  satisfying  the  following  conditions 

ceC^  (M^)  ,c{x)  =  lioT  X  e  (126) 

c{x)  >  Co  =  const.  >  0,\fx  E  G.  (127) 

In  addition,  we  assume  that  there  exists  a  point  Xq  G  D  such  that 

[x  —  xo,Vc~‘^  (x))  >  0,\/x  E  G.  (128) 

Condition  (128)  can  be  interpreted  as  the  monotonicity  condition  of  the  function  {x) 
along  rays  of  straight  lines  passing  through  the  point  {xq}  ,  and  {xq}  is  the  beginning  of 
any  such  ray.  We  use  (128)  quite  essentially.  First,  we  use  it  in  Lemma  1.  This  lemma  is 
applied,  in  turn  to  establish  the  exponential  decay  with  respect  to  the  time  t  of  the  solution 
of  the  Cauchy  problem  for  the  acoustic  equation  in  the  time  domain  using  results  of  the 
book  [68].  Second,  (128)  guarantees  the  existence  of  the  Carleman  estimate  for  the  operator 
—  (?  (x)  A,  see  Theorem  2.6  in  [13].  This  Carleman  estimate,  in  turn  allows  us  to  apply 
Theorem  3.1  of  [13]  in  the  hnal  step  of  the  proof  of  Theorems  9.5  and  9.6  below. 

Lemma.  Assume  that  conditions  (126)-(128)  are  in  place.  Then  the  family  of  geodesic 
lines  generated  by  the  function  c(x)  satisfies  the  non-trapping  property  in 
Consider  the  function  g  {x)  satisfying  the  following  conditions 

qEG'^  (M^)  ,g{x)  =  0  in 

g{x)  ^0,xE  S. 

Ag  (x)  7^  0,  Vx  G  Gi. 

Let  the  number  /c  G  M.  Consider  the  following  problem 

k‘2 

Au  +  „  ,  =  -g  (x)  ,x  E 

[x) 

3 

7  -^dxu(x,k)—iku(x,k)  =  o(l)Ax\^oo.  (133) 

\x\ 

j=i 

Equation  (132)  is  the  acoustic  equation  in  the  frequency  domain.  Just  as  in  the  previous 
section  5.3.2,  we  now  refer  to  Theorem  6  of  Chapter  9  of  the  book  [68],  Theorem  3.3  of  the 
paper  [67]  as  well  as  to  Theorem  6.17  of  the  book  [45].  Combining  these  results  with  Lemma, 


(129) 

(130) 

(131) 

(132) 
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we  obtain  that  for  each  /c  G  M  there  exists  unique  solution  u  {x,  k)  G  (7®+“  (M^) ,  Va  G  (0, 1) 
of  the  problem  (132),  (133). 

Phaseless  Inverse  Problem  5  (PIP5).  Suppose  that  the  function  c{x)  satisfying 
conditions  (126)-(128)  is  unknown  for  x  E  Q  and  known  for  x  G  Assume  that  the 

following  function  /s  {x,  k)  is  known 

/s  {x,  k)  =  \u  {x,  /c)  I  ,  Vx  G  S',  V/c  G  (a,  b) .  (134) 

Determine  the  function  c  {x)  for  x  E  kl. 

Theorem  9.5.  Consider  PIPS.  Let  the  function  g{x)  satisfies  conditions  (129)-(131). 
Consider  two  functions  Ci  {x) ,  C2  {x)  satisfying  conditions  (126)-(128)  and  such  that  Ci  {x)  = 
C2  (x)  =  c(x)  for  X  E  M^\D.  For  j  =  1,2  let  Uj{x,k)  E  (M^) ,  Vo  G  (0,1)  be  the 

solution  of  the  problem  (132),  (133)  with  c(x)  =  cj  (x).  Assume  that 

\ui  {x,  k) I  =  \u2  {x,  k)\  ,\/x  E  S,Wk  E  (a,  b) .  (135) 


Then  Ci  {x)  =  C2  (x) . 

IPl  is  about  the  case  when  the  modulus  of  the  total  wave  held  is  measured  for  x  E  S,  k  E 
(a,  b).  Consider  the  function  uq  (x,  k) , 

,  , ,  f  exp  {ik  \x  -  ^1)  _ 

^o{x,k)  =  j  giOdt 

G 

This  function  is  the  solution  of  the  problem  (132),  (133)  for  the  case  c  {x)  =  1.  Since  c{x)  =  1 
for  X  E  M^\G,  then  one  can  consider  uq  {x,  k)  as  the  solution  of  the  problem  (132),  (133)  for 
the  background  medium.  Hence,  the  function  Ug  [x,  k)  =  u  {x,  k)  —Uq  [x,  k)  can  be  considered 
as  the  wave,  which  is  scattered  due  to  the  inhomogeneous  structure  of  the  coefficient  c  [x) 
for  X  E  G.  This  is  our  motivation  for  posing  Phaseless  Inverse  Problem  6. 

Phaseless  Inverse  Problem  6  (PIP6).  Suppose  that  the  function  c{x)  satisfying 
conditions  (126)-(128)  is  unknown  for  x  E  Ll  and  known  for  x  E  Let  Us{x,k)  = 

u  {x,  k)  —  Uq  {x,  k) .  Assume  that  the  following  function  /g  (x,  k)  is  known 

/e  (x,  k)  =  1^5  (x,  k)  \  ,'ix  E  S^'ik  E  (a,  b) . 

Determine  the  function  c  {x)  for  x  Ekl. 

Theorem  9.6.  Consider  PIP6.  Assume  that 

C  (x)  7^  1,  Vx  G  S. 

Let  all  conditions  of  Theorem  5.5  hold,  except  that  (130)  is  not  imposed.  In  addition,  let 
(135)  be  replaced  with 


{x,  k) I  =  \us^2  {.X,  k)\  ,'ix  E  S,'ik  E  (a,  b) , 
where  Ugj  {x,  k)  =  uj  {x,  k)  —  uq  {x,  k)  ,j  =  1,  2.  Then  ci  {x)  =  C2  {x) . 
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10  Presentations 


Below  is  the  list  of  our  presentations  of  above  results  during  the  reportage  period. 

A.  Presentations  to  engineers  of  the  Army  Research  Laboratory  (ARL): 

1.  Presentation  at  the  Army  Research  Laboratory  (ARL)  to  Mr.  Brian  Burns,  an  engineer 
of  ARL,  Dr.  Joseph  D.  Myers,  Program  Manager  of  Numerical  Analysis  Program  of 
ARO  and  to  Drs.  Anders  Sullivan  and  Lam  Nguyen,  engineers  of  RF  Signal  Processing 
and  Modeling  Branch  of  Sensors  and  Electron  Devices  Directorate  of  ARL  and  to 
(2012). 

2.  Presentation  to  Mr.  Brian  Burns  and  to  Dr.  Joseph  D.  Myers  during  their  visit  of  our 
campus  (2012). 

3.  Presentation  to  Drs.  Joseph  D.  Myers  and  Anders  Sullivan  during  their  visit  of  our 
campus  (2012). 

4.  Presentation  over  the  phone  to  Mr.  Brian  Burns  and  Drs.  Joseph  D.  Myers  and  Anders 
Sullivan  (2013). 

B.  Other  presentations: 

1.  Colloquium  talk  at  Department  of  Mathematical  Sciences  of  University  of  Delaware. 

2.  Colloquium  talk  at  Department  of  Mathematics  of  University  of  Virginia. 

3.  Colloquium  talk  at  Department  of  Mathematics  of  Georgia  Institute  of  Technology. 

4.  Colloquium  talk  at  Department  of  Mathematics  of  Wichita  State  University  (Wichita, 
Kansas). 

5.  Inverse  Problems  workshop  in  Ecole  Polytechnique,  Paris. 

6.  Plenary  speaker  at  the  conference:  “Inverse  Problems  Problems:  Modelling  and  Sim¬ 
ulation”,  Turkey. 

7.  Plenary  speaker  at  the  conference  on  Inverse  Problems,  Albi,  France. 

8.  Conference  “Applied  Inverse  Problems”,  Korea,  organizer  of  a  minisimposium  and 
presenter  of  a  talk. 

9.  Presentation  at  the  Conference  dedicated  to  the  30st  anniversary  of  the  journal  “Inverse 
Problems”,  Bristol,  UK. 

10.  Presentation  at  the  6th  International  Conference  on  Advanced  Computational  Methods 
in  Engineering,  Gent,  Belgium. 

11.  Plenary  speaker  at  International  Conference  on  Inverse  Problems  and  Optimal  Control, 
Hong  Kong. 
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11  Book  and  Survey  Papers 

The  book  [1]  is  published  by  Springer,  the  World  most  prestigious  publisher  of  the  scientihc 
literature.  Results  of  Beilina  and  Klibanov  since  their  first  joint  publication  [25]  in  2008  and 
up  to  2012  are  summarized  there.  The  book  discusses  both  the  theory  and  numerical  studies 
of  the  globally  convergent  numerical  method  of  section  3  and  of  the  adaptivity  of  section  4. 
Although  this  book  was  published  only  two  years  ago,  it  was  cited  more  than  97  times  by 
now,  see 

http:/ /scholar. google. com/citations?user=pFmp7LMAAAAJ&hl=en. 

The  paper  [7]  is  the  survey  of  the  relaxation  property  for  the  adaptivity  technique  (section 
4)  .  The  paper  [13]  is  the  survey  on  the  applications  of  Carleman  estimates  to  inverse 
problems,  starting  from  the  work  [37]  in  1981.  It  is  worthy  to  point  out  that  both  globally 
convergent  numerical  methods  for  CIPs,  presented  above,  have  roots  in  the  idea  of  [37]. 
Currently  there  are  more  than  50  publications  of  authors,  other  than  the  ones  of  [37],  which 
cite  this  paper.  In  accordance  with 

http:/ /scholar. google. com/citations?user=pFmp7LMAAAAJ&hl=en, 
the  paper  [37]  was  cited  more  than  314  times. 


12  Publications 

Publications  [l]-[23]  acknowledge  the  support  of  this  grant.  The  paper  [27]  is  a  featured 
article  of  Inverse  Problems. 
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